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ABSTRACT 



Transient difiuser flow in an HF/DF laser is analyzed via a simplified model. A radial laser 
configuration is assumed where the diffuser is idealized as a ring- symmetric channel of uniform width. 
The outward radial flow following an abrupt start-up of diffuser inflow is computed numerically using 
a GRP-type code for integrating the Euler equations, which include an ongoing reaction of hydrogen 
recombination. A new kinetic scattering effect is suggested, resulting from the recoil of an HF or DF 
molecule assisting as third body in the hydrogen recombination reaction. The main conclusion is that 
transient diffuser plume flow effects do not drastically aggravate the potential for contaminating 
backflow (HF + DF) caused by either thermal or kinetic backscattering. 



11 



ACKNOWLEDGEMENT 



This work was conducted as part of a laser exhaust study under the cognizance of Distinguished 
Professor Allen E. Fuhs. 1 deeply appreciate and wish to thank Professor Fuhs for his continuous 
support and guidance. 



TABLE OF CONTENTS 



1. INTRODLCTION 1 

2. THE DIFFUSER FLOW MODEL 4 

2.1 Radial Laser Configuration 4 

2.2 Typical Laser Flow 4 

3. THE GASDYNAMIC FLOW MODEL 7 

3.1 The Governing Equations 7 

3.2 Hydrogen Recombination 8 

4. RESULTS AND DISCUSSION 10 

4.1 Thermal Backscattering 10 

4.2 Kinetic Backscattering 13 

4.3 Concluding Remarks 15 

5. REFERENCES 18 

APPENDIX A. The HFL Code 20 

A.l Features Specific to DifTuser Flow Model 20 

A. 2 Code Listing 22 

6. DISTRIBUTION LIST 43 



LIST OF FIGURES 



Figure l-l. Ring-Symmetric HF/DF Laser Exhaust Plume 3 

Figure 2-1. Radial Configuration of HF DF Open-Loop Laser 6 

Figure 4-1. Mach Number at Diffuser Exit for Reactive Transient Flow 16 

Figure 4-2. Rotation Angle for Non-Steady Plume (Overestimate) 17 



IV 



NOMENCLATURE followed by units (if any) 



A cross-section area (of diffuser or nozzle) 

A^, Avogadro's number 6.022 X 10"^ (molecules kmole) 

C sound speed (m sec'*) 

E total energy per unit volume (MJ m'^) 

fy mole fraction of species X at diffuser entrance 

gp rate of hydrogen recombination collisions (with HF or DF) per unit volume (m*^ sec'*) 

hydrogen recombination rate constant [(mole cnr^)'" sec'*] 
k| modified hydrogen recombination rate constant [(kg nr^)'" see'*] 

M Mach number 

n number density (molecules m^) 

P pressure (Pa) 

formation energy of H, (FJ mole) 

Qd£t hydrogen recombination energy release parameter (MJ/kg) 
t time (ms) 

T temperature (K) 

U flow velocity (m sec) 

W average molecular weight (kg kmole) 

X radial coordinate (in diffuser) 

Z normalized mole fraction of H (Z= 1 at diffuser entrance) 

a turning angle of flow velocity in a corner expansion to vacuum 

Aa turning angle increment for transient corner expansion 
y specific-heat ratio 

p density (kg m'^) 

0 inclination of flow velocity vector (to x-axis) 

r the fraction (y+ l)/(y-l) 

[X] molar density of species X (kmole m'^) 
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stagnation conditions 
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nozzle (difiiiser) exit conditions 


02 


dilTuser entr\’ conditions 


()t 


Eulerian time derivative 


Ox 


Eulerian space derivative 


• 

0 


Lagrangian time derivative 



INTRODUCTION 



In recent years a study of the gasdynaniics of space-based IIP DF laser exhaust plume was 
conducted at the NPS. The aim of this program was to identify phenomena that give rise to a 
potentially contaminating backflow of corrosive species (UP. DF, F), and to estimate the magnitude 
of the flux arriving at the spacecraft (I, 2, 3, 4] . The laser was envisioned as a cylindrical spacecraft 
having a centrally located nozzle of ring-symmetry. This is an idealization of a more general zero- 
thrust exhaust configuration for an open loop IlF/DF laser (Fig. 1-1). 

In former studies (2, 3] we focused on the contribution of thermal self-scattering to a 
contaminating backscattered flux of corrosive molecules (HF, DF) from the exhaust plume. It was 
shown that this flux emanates primarily from the lip-centered rarefaction fans that flank the e.xhaust 
plume, and it was estimated that at a sufficiently high exit .Mach number (e.g. Mj = 4), the thermally 
backscattered flux of HF + DF is utterly negligible. 

The purpose of this report is to present a study of the effect of diffuser start-up flow, including the 
ongoing recombination of atomic hydrogen, on the thermally backscattered flux. Our main 
conclusion is that by designing for a sufficiently high steady exit Mach number (e.g. Mj = 4), the level 
of thermally backscattered flux of corrosive molecules (HF-i-DF) would be negligible even when the 
combined effects of hydrogen recombination and transient flow are considered. 

The diffuser is simplified as a cylindrically expanding channel of uniform width, and the flow is 
idealized as inviscid expansion of a perfect gas with an ongoing reaction of hydrogen recombination. 
The flow is computed as one-dimensional time-dependent. 

A fully 2-D time-dependent computation of an emerging exhaust plume is outside the scope of the 
present laser exhaust study. Rather, we resort to semi-quantitative arguments that lead to an upper- 
bound estimate of the effect of transient flow on thermal backscattering, by showing that the 
transient turning angle in an expansive corner flow can be estimated as higher than the corresponding 
steady flow turning angle. It is subsequently suggested that an overestimate to the the backscattered 
flux from a transient plume is obtained by considering it as a steady plume whose exit velocity vector 
is pre-rotated so that its limiting (vacuum) streamline coincides with that defined by the transient 
turning angle. 
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The plan of this report is the following. In Ch. 2 we present our radial (1-D) dilTuser model based 
on the linear configuration of the TRW test laser (5] , and we propose a "typical case" of HF DF 
laser flow based on one of those tests. Ch. 3 is devoted to the gasdynamic governing equations, 
including the hydrogen recombination rate. The resulting diffuser model was implemented in a 1-D 
Euler code (named IIFL) which utilizes the GRP scheme for integrating the conservation laws of 
compressible flow [6] . The HFL code is given in Appendix A. The results of a typical diffuser start- 
up flow are presented in Ch. 4, along with a discussion and analysis of the effects of transient flow' 
and uncertainty in the hydrogen recombination rate on thermally backscattered flux of HF+ DF. We 
also consider kinetic backscattering of HF or DF molecules caused by third-body recoil in the 
hydrogen recombination reaction. It is shown that the combined contribution of these effects to 
backscattered flux remains negligible in the typical case. Chapter 4 ends with some concluding 
remarks, and references are brought in Ch. 5. 
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SPACE - BASED HF LASER 



Figure 1-1. 
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Ring-Symmetric HF DF Laser Exhaust Plume 
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RING - SYMMETRIC PLUME 



2. THE DIFFUSER FLOW MODEL 



Our difTuser model is based on a schematic radial laser configuration as shown in Fig. 2-1. The 
major components in this design concept are those present in an experimental HF, DF open loop laser 
tested at TRW [5] . In this chapter we describe the radial configuration including the difTuser 
(Section 2.1), then we present a typical HF DF laser flow based on one of those TRW tests 
(Section 2.2). 



2,1 Radial Laser Configuration 



We assume that an open loop HF/DF laser of the type tested at TRW [5] , can be rearranged in a 
radial (rather than linear) sequence, where the flow begins at the hub and proceeds in an outward 
radial direction. Referring to Fig. 2-1. the major components of this configuration are : 

(a) Combustion Chamber : Deuterium and fluorine burn with excess fluorine, resulting in hot 
gaseous mixture where fluorine is virtually completely dissociated. 

(b) Nozzle Cascade ; Rapid expansion to supersonic flow leaves atomic fluorine concentration 
effectively frozen. 

(c) H, + He Injection: Mixture is injected between the supersonic streams of DF + F emerging 
from the cascade. 

(d) Vlixing and Lasing : The lasing is from vibrationally excited HF molecules produced by direct 
reaction between H, and F. As a by-product, one H atom is produced for every HF molecule. 

(e) DifTuser Entrance : This point marks the end of the lasing process. From this point on the 
flow is just an exhaust to be discarded safely. 

(0 Radial DifTuser : The purpose of the diffuser is to raise the flow Mach number at the exit so 
that no appreciable backflow from the exhaust plume will take place [2, 3, 4] . It should be 
noted that a desirable diffuser area ratio can be achieved at lower radius ratios by letting the 
diffuser expand in the axial direction. 

2,2 Typical Laser Flow 

We chose as a typical laser flow one of the tests conducted in the TRW experiments [5] it is test 
III in Table 5 (p. 91) of that reference. The mole fraction and corresponding average molecular 
weight and specific-heat ratio (W and y) at the diffuser entrance for this test w'ere : 
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Mole fractions 



Average molecular weight 



fn = 


.091 


*^IIF 


.091 


f|i, - 


^DF ^ 


.135 


*^11 e ^ 


.579 





W = 7.27 



Specific-heat ratio 



y = 1.5-4 



( 2 . 2 - 1 ) 



The additional data required for the computational model is the flow variables at the difluser 
entrance. Some of these variables have been measured directly or indirectly and the remaining 
variables can be evaluated from standard isentropic flow relations for compressible flow of an ideal 
gas. The data given in the TRW tests report [5] are : 



Difluser entrance cross-section is : 1" x 7" 



Mass flow rate is : 14.99 (gr sec*’) 

Stagnation temperature is ; 1400 (K) 

Flow velocity is ; 2300 (m sec'*) 



( 2 . 2 - 2 ) 



Flow density is now obtained as the ratio of the specific mass flow rate and the velocity. Then 
Mach number is extracted from the standard relations between Mach number, velocity and stagnation 
temperature. The values of these variables are : 



M, = 2.26 

p2 = 1.44 X 10'^ (kg m'^) (2.2-3) 

P 2 = 9.72 X 10-^ (MPa) 

The flow at the diffuser entrance is thus fully specified for the selected typical case. In the next 
chapter we take up the matter of computing the transient expansion following an abrupt start-up of 
the diffuser inflow. 
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Figure 2-1. Radial Configuration of HF'DF Open-Loop Laser 
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3. THE GASDYNAMIC FLOW MODEL 



The transient expansion of the lasing products following an abrupt start-up of inllow at the 
difluser entrance, is idealized as an inviscid compressible flow of a mixture of perfect gases with one 
ongoing chemical reaction — that of hydrogen recombination. The computational model thus calls 
for a formulation of the governing equations and for a numerical scheme capable of integrating them 
in time and space. 

In this chapter we describe the governing equations for the diffuser flow (Section 3.1), following by 
a simple approximation chosen for the hydrogen recombination rate as a three body reaction (Section 
3.2). The numerical scheme employed for the solution is of the Generalized Riemann Problem (GRP) 
type [6] . This scheme has been implemented in a code which was adapted to the diffuser How case. 
The code (named HFL) and some brief description of features introduced to treat hydrogen 
recombination are given in Appendix A. In a recent report (7) a very similar version of this code 
(without chemical reaction) was described in considerable detail. 



3.1 The Governing Equations 

The expansion of lasing products through the diffuser is governed by the Euler equations for a 
gaseous mixture of ideal gases with ongoing hydrogen recombination reaction. We write these 
equations in the quasi-one-dimensional format of flow in a stream tube of varying cross-section area 
A(.\), but in actual computations we set A(x) = .\, thereby reducing the equations to the cylindrical 
case. (The code HFL, however, can accept any smooth A(x)). 

We chose to simplify the designation of atomic hydrogen mole fraction by normalizing it as 
fHZ(x.t), where fj_j is the mole fraction at the diffuser entrance (it is constant) and Z(x.t) denotes the 
degree of dissociation (Z= 1 is maximum concentration of H at the diffuser entrance, Z = 0 is 
complete recombination). As an approximation, we assume that W and y are constant throughout. 
In the typical case (2.2-1), y and W can vary by as much as 1% and 4 % respectively at full 
recombination. Neglecting that variation is consistent with the preliminarx' nature of the present 
diffuser flow analysis. 

The governing equations are the three standard conservation laws (mass, momentum and energy), 
and an additional species conservation law for atomic hydrogen. The conservation laws are 
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augmented by an equation of state (ideal gas) and by a rate law for hydrogen recombination. The 
full system of equations is : 



Mass 


lAp], + [ApUj^ =0 


Species H 


[ApZ], + [ApUZl^ = ApZ 


•Momentum 


[ApU]^ + [pV\ AP^ =0 


Energy- 


[AE], + [A(E + P)U]^ =0 


Equation of State 


P(p,E,Z,) = (y-l)(E - pU’‘/2 - (pf„Z/W)Q„l 


Recombination Rate 


Z = Z(p,P,Z) 



where A is a function solely of x, and p, U, P, E, Z are all functions of (x.t). The normalized rate 
function is Z. Its derivation and explicit expression are given in Section 3.2. 

We note that the energy equation contains the heat released by hydrogen recombination (Q^ per 
mole of H) in an implicit way, by defining the total energy as including the latent heat due to 
recombination, in addition to internal and kinetic energy terms. 



3.2 Hydrogen Recombination 

The hydrogen recombination is commonly assumed to be a three-body reaction [8] , with the 
following rate law ; 

H + H + M H 2 + M 

IH 2 I = (3.2-1) 

= 7 X 10 ’^ l(mole;cm^)‘“ sec"'] 
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where by M we denote a third molecule, which in our case may be any molecule in the gaseous 
mixture (including H). 



The value we chose for is an upper limit of a range of values recommended by Kondratiev [h] 
for the reaction H + H + He (from 2 x 10*^ to 7 x lO*^). since in the typical case helium 
constitutes about 58% of the molecules in the lasing products. 

The data compiled by Kondratiev [9] indicates a comparable range of variation with temperature 
as with the third species M. We later introduce uncertainty factors of 10 and 100, which in all 
likelihood more than reflect the uncertainty in the reaction rate for hydrogen recombination. 



It is convenient to redefine the reaction rate in terms of flow density p and normalized H 
concentration Z. The modified rate constant kj is given by : 

Z = - k, p- Z“ 

(3.2-2) 

k, = 2kof„/w2 



where the factor 2 is due to the fact that the rate of depletion of H is twice the rate of production of 
H,. In deriving (3.2-2) we made use of the relations = IH],/1M], and [M]=p/W, where index 2 
denotes the dilTuser entrance. 



The heat of formation released per mole of recombined H 2 is 435.783 (kJ mole) according to [10] 
(page F-179 in that reference). In the code we use the energy parameter per unit mass of the 
gaseous mixture Qq£j, defined as : 



Qdet Qh 

= 435.783 / 2 (KJ;mole) 



(3.2-3) 
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4. RESULTS AND DISCUSSION 



Several numerical computations of the transient dilTuser (low were performed using the code HFL 
(Appendix A). The flow was started by an abrupt inflow at the diffuser inlet; the radial diffuser 
model and the typical HF DF laser case (Sections 2.1 and 2.2) were assumed. Due to uncertainty in 
hydrogen recombination rate, three rate levels were assumed : the nominal level (3.2-1), a tenfold 
increased rate and a hundredfold increased rate. 

We are primarily concerned with the thermally backscattered flux arriving at the surface from the 
exhaust plume. Our goal is to show that by choosing a diffuser with a sufhciently high steady exit 
Mach number, the backscattered flux can be negligibly low even when the combined elTects of 
hydrogen recombination and transient exit flow are taken into account. Put in other words, we 
present a simplified analysis of the transient exhaust flow, which enables the establishment of a 
reasonable margin on a frozen/steady design for low backscattered flux. This analysis is presented in 
Section 4.1 below. 

The continuation of the hydrogen recombination reaction into the plume may give rise to a 
contaminating backflow of a different kind. The velocity imparted to the "third body" molecule in the 
hydrogen recombination process (Eq. 3.2-1) may be large enough to overcome the radial flow velocity 
component, resulting in molecular backflow. This effect is discussed in Section 4.2, and it is shown to 
be potentially capable to give rise to an HF/DF deposition rate of the order of 1 molecular 
monolayer per hour. 



4.1 Thermal Backscattering 

Since in steady flow thermally backscattered flux is linked primarily to exit Mach number, we focus 
our attention on the time-history of the exit .Mach number obtained from the diffuser start-up 
computations mentioned above. Two major features are noted in these results (Fig. 4-1). 

The first feature is the monotonic decrease in exit Mach number with time, as the flow within the 
diffuser undergoes transition from an initial "cloud expansion" mode to a steady expansion flow in a 
channel of increasing cross-section area. A steady exit flow seems to be established after about 
2 (ms) from the start-up instant. The second feature is related to the heat released in the flow by 
hydrogen recombination. As the recombination rate is increased, so does the time-asymptotic value 
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of the exit recombination fraction. The recombination fractions in the three cases (Fig. 4-1) are I‘*o. 
7% and 46%; the corresponding exit Mach numbers are 4.0, 3.69 and 2.81. This trend is in 
qualitative agreement with the well known phenomenon of decrease in .Mach number as a result of 
heat addition to a steady supersonic flow in a channel of uniform cross-section [11] . 

Estimates of recombination rates are notoriously uncertain. How can a reasonable uncertainty 
factor be established? We do not know of an estimate of this factor for the How of HF DF laser 
products. However, in another case of flow involving hydrogen recombination - that of hydrazine 
rocket motors - a recent study [12] recommends a hydrogen recombination rate of 2.8 x lO’^/T (in 
c.g.s. K units) and an uncertainty factor of 30. Assuming the relatively low temperature ofT=300 
(K) and multiplying by 30, we get a recombination rate of 3 x lO'^, versus ?x 10*^ in our tenfold 
case. This demonstrates that the tenfold case is already reasonably high; just the same, we also retain 
the hundredfold case in the upcoining discussion. 

The exit Mach number time-histories (Fig. 4-1) demonstrate that the nominal diffuser flow is 
nearly frozen (exit recombination fraction 1% and Mach number 4); even in the tenfold case the flow 
is not dominated by hydrogen recombination (exit recombination fraction 7% and .Mach number 
3.69). This conclusion is in agreement with findings of the TRW HF/DF laser study [5] . It takes the 
unrealistically high hundredfold increase in hydrogen recombination rate to produce a decisive change 
in the steady diffuser flow (exit recombination fraction 46% and .Vlach number 2.81). 

One role of the diffuser is to e.xpand the laser exhaust to a sufficiently high exit Mach number, so 
that the thermally backscattered flux of corrosive molecules (HF, DF) is negligibly small. Let us 
consider the nominal case (Mj = 4) and denote by "reference flux" the combined HF+DF 
backscattered flux arriving at a point located at 0.1 (m) from the nozzle lip (Fig. 1-1). Using the code 
RINGED based on the breakdown surface model [2] , the nominal reference flux was evaluated as 
3.1 X 10^ (molecules m'" sec"'). This flux level corresponds to a surface deposition rate (assuming the 
sticking coefficient equals 1} of about 10"'® molecular monolayers per hour, which is utterly negligible 
(a reasonable estimate of total operating time would not exceed several hours). 

Observing that a steady exit flow' has been established in all three cases by about t = 2 (ms) (see 
Fig. 4-1), we suggest that an approximate estimate of the backscattered flux would be obtained by 
assuming a steady flow at the diffuser exit and the exhaust plume. The RINGED computations of 
the increase hydrogen recombination rate cases yielded the following results. In the tenfold case 
(Mj = 3.69) the flux was 10' times larger than the nominal reference (lux. In the hundredfold case 
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(M| = 2.81) the flux increase was by a factor of 10^. Even in the unrealistically high hundredfold case, 
that flux level corresponds to about 10'^ monolayers per hour, which is still negligible. 

What other uncertainties could affect the foregoing conclusion? First we consider the effect of 
variations in stagnation temperature and density relative to the nominal case, then we take up the 
matter of transient plume flow. 

E.xit stagnation temperature was observed to vary' in the HFL diffuser start-up computations by a 
factor of no more than about 1.3. Since all velocities in the breakdown surface model are normalized 
by speed of sound [2] , the flux is proportional to the square root of the stagnation temperature. 
Also, it has been shown (2] that a change in stagnation density can produce an inversely proportional 
change in flux. Stagnation density was reduced by a factor of about 1.3 in the tenfold case and about 
2.S in the hundredfold case. In either case the potential effect on flux is small relative to changes 
brought upon by reduced exit Mach number, so we attribute the drastic change in flux following an 
increased hydrogen recombination rate primarily to the reduced exit Mach number. 

We now address the effect of transient plume flow. Since the exit Mach number is monotonically 
decreasing with time (Fig. 4-1), and since other flow variables do not appreciably affect the flux from 
a steady plume, we suggest that an upper bound on the effect of transient plume flow can be 
established by linking it to a "transient turning angle", larger than that corresponding to a steady flow 
with the same exit Mach number. Obviously, if a Prandtl-Meyer flow pattern is derived from a flow 
which e.xits the nozzle at an angle lower than 90", its limiting turning angle would bring it closer to 
the spacecraft surface, and with a lower outward radial velocity to overcome, more molecules would 
be thermally backscattered. Thus, we regard the transient flow around the nozzle lip as represented 
by an equal exit Mach number steady flow, with some initial exit turning angle. 

The justification for this model is the following "upper bound" estimate for the transient turning 
angle of an emerging exhaust plume. The exit velocity vector Uj is visualized as rotating by a gradual 
turning of the plume/vacuum interface until it reaches the normal velocity (2/(y-l))C| (Fig. 4-2). 
which corresponds to a plane-wave expansion into vacuum. The corresponding turning angle is 
a = (2/(y-l))C|/U| = (2/(y-l))/Mj. It is an overestimate of the turning angle since in a steady corner 
expansion (and also in a nonsteady one) the magnitude of the velocity vector increases as it rotates, 
in order to conserve stagnation enthalpy (or even increase it in nonsteady flow about an expansive 
corner). 
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It is of interest to note that ct as given above is identical with the first term of a power series 
expansion in of the standard Prandtl- Meyer expression for the limiting turning angle (13| . Let 

us denote by A« the increment of the unsteady turning angle relative to the corresponding steady one 
(Prandtl-Meyer). This (positive) increment is given by : 

Aa = (2/(y-l))/M| - | F ^ arctan[r ^ (Mj'-i) ^ ] — arctan[(M,‘ — l) ^ } 

( 4 . 1 - 1 ) 

r = (y+i)/(y-i) 



We now argue that the flux ratio between a steady flow and the same flow with exit velocity 
rotated by Aa is an upper-bound estimate to the effect of transient plume flow with the same exit 
variables. We note that by pre-rotating the steady corner How. we make the limiting streamline in the 
steady case, coincide with the nonsteady plume'vacuum interface deemed to have rotated through the 
angle a. 

The stage is now set to estimate the transient flux increase factors in the three hydrogen 
recombination rate cases (exit .Mach numbers 4, 3.69 and 2.81). Using Eq.(4-1) we get the turning 
angle increments of Aa= 4.1”, 5.1" and 10.6" respectively. Re-computing the reference flux (code 
RINGED) in these cases with initial exit angle of 90"- Aa, we get the fiux increase factors of 2.7, 3.5 
and 10 respectively. Even in the worst case (hundredfold), the flux would increase from 10'^ to lO'" 
monolayers per hour, which is still negligibly small. 

We also notice that the fiux increase factor resulting from the decrease in exit .Mach number due 
to higher hydrogen recombination rate, is much larger than the factor representing the effect of 
transient plume flow in the same case. Thus, in the tenfold case these factors are 10" versus 3.5, and 
in the hundredfold case the ratio is even higher ; 10^ versus 10. We conclude that the transient effect 
is quantitatively secondary to the effect of exit .Mach number, which is thus established as the major 
parameter for designing an exhaust with a negligible level of thermally backscattered HF+ DF flux. 
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4.2 Kinetic Backscattering 



The role of the third body in the hydrogen recombination reaction (3.2-1) is to absorb the energy 
released by the recombination process while maintaining the combined pre-collision momentum of 
the participating molecules. Assume the momentum added to [MJ is q, then the momentum added to 
H, is — q. Denote by itIj the mass of H,, by ni^ the mass of HF and by e the energy released due to 
the recombination. Then q is determined from the following energy conservation relation ; 

q-|2m^ + qV^m, = e 

(4.2-1) 

ntj = 2 m 2 = 20 e = 4.36 x 10* (Joules kmole) 

U = q/m, = 2000 (m, sec) 

The values used in the computation above are per kmole, but the result is the same as using values 
pertaining to single molecules. Since the velocity obtained upon expansion to zero-pressure in the 
typical case is about 3000 (m sec), our results imply that a turning angle in excess of arccos(2,'3) = 48° 
is needed in order to enable some kinetically backscattered molecules to reach the spacecraft. 

In the nominal case the turning angle is 49" and in the tenfold case it is 52". Also, some 
kinetically backscattered molecules can originate from within the rarefaction fan. since the exit 
velocity in the typical case is only about 1500 (m sec'*) (versus 2000 (m sec’’) velocity increment to 
kinetically scattered HF molecules). Consequently, this effect may contribute to spacecraft 
contamination, and its magnitude should be estimated. It will become progressively more significant 
at lower exit Mach number, and thus may be an additional factor in determining a minimal value of 
this flow variable. 

While we do not presently attempt at an exact integration of the flux of kinetically backscattered 
molecules arriving at each surface point, we propose the following overestimate of its magnitude. 
Consider the outgoing flux of kinetically scattered HF molecules from a half-space of quiescent 
uniform gas having the thermodynamic properties of the nozzle exit in the typical case. Consider the 
hydrogen recombination rate equation (3.2-1), and substitute pi] =fj^n'A^. This equation then reads 
as the volume rate of hydrogen recombinations (denoted as below). In a slab of quiescent gas. half 
the kinetically scattered molecules have an outward pointing velocity vector in one direction. The 
probability of collisionless passage out of the slab is exp( — x/X.), so that the outgoing flux F is given 
by : 
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F = g, X/2 

(4.2-2) 

"c ^ *^0 ^Dp) 

Using typical case exit conditions nj = 2.81 x 10"^ (nU'^) and mean free path Xj = 1.28 x lO"* (m). 
we get F = 5.14x lo'^’ (m'" sec''), or about 12 monolayers per hour. Since only a small fraction of 
this flux actually arrives at spacecraft surface points, the arriving flux is in all likelihood well below 1 
monolayer per hour. If this flux level is deemed significant, an accurate integration of kinetically 
scattered flux emanating from the exhaust plume should be performed. The details of such scheme 
are analogous to the first-collision ambient scattering model (4J , since both elTects result in a source- 
like distribution throughout the plume, and both involve a steric factor (one half in the quiescent gas 
above) related to the vector addition of flow velocity and velocity imparted through ambient collision 
or assistance in hydrogen recombination. Therefore, the two effects can conveniently be unified 
under a common framework. 



4.3 Concluding Remarks 

From the foregoing analysis and discussion we conclude that a design for low (negligible) level of 
thermally backscattered flux of HF + DF can be accomplished by assuming steady frozen exhaust 
flow. The effects of transient diffuser plume flow and hydrogen recombination may typically increase 
that flux by a factor of 10". Given the sensitivity of flux to exit .Mach number, an adequate margin 
can readily be achieved by raising the exit Mach number. In the typical case presented here, an exit 
Mach number of 4 has been shown to be adequate. 

It should be noted that transient effects originating from points upstream of the diffuser entrance 
were not considered here, as they depend on details of the system construction and operating 
sequence. 

The kinetic scattering effect has been pointed out as an additional potential source of 
contaminating backflow. It is the result of hydrogen recombination collisions assisted by an IIF or 
DF molecule as third body. The HF/DF molecule receives part of the kinetic energy release by the 
hydrogen recombination. 
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Figure 4-1. .Mach Number at DifTuser Exit for Reactive Transient Flow 
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APPENDIX A. The HFL Code 



A.l Features Specific to Diffuser Flow Model 

The purpose of this appendix is to provide the listing of the HFL code used in the diffuser flow 
computations. An almost identical code vv'as recently described in some detail [7] , the major 
modification for HFL being the addition of hydrogen recombination. 

The highlights of the HFL version to the GRP code [71 are as follows (numbers in parenthesis 
refer to statement numbers which are prefi.xed by HFL in the listing). 

(a) Data (NETUN.M) : 

The data which relates specifically to the diffuser flow with hydrogen recombination is given in 
statements (194 - 239). Some parameters are stored in /DIFFUS, (129) for use in BEGIN and 
SAFAE. 

(b) Initial Conditions (BEGIN) : 

The initial conditions are ideally an empty diffuser. As an approximation we set small pressure 
and density (270 - 271) and a positive velocity (277) as initial values. The velocity reduces the 
intensity of the shock driven into the dilute initial gas by the abrupt start-up of inflow at the 
diffuser entrance (1 = 2). Note that the initial value of the total energy E(l) is augmented by the 
heat of formation of hydrogen (306). The heat term is subsequently subtracted from the total 
energy in order to compute the pressure (see (d) below). 

The grid is defined as follows. The number of cells (L-2) is L-2 = 50 (24). The diffuser e.xtends 
from X0 = 0.625 to XI = 3.125 (2S3 - 2S4), so that DX = 0.05. Since the external radius of the 
spacecraft is just 2.5, this grid consists of 40 cells within the diffuser and 10 cells in an extended 
segment. The purpose of this extended segment is to minimize any propagation of error due to 
imperfection of the outflow boundary condition (see (c) below). Since the ambient pressure is 
zero, this additional segment does not affect the flow at X 2.5. Consequently, the results for 
the diffuser exit flow were read from cell 1 = 39 whose midpoint is X = 2.5. 

(c) Boundary Conditions (SAFAE) : 

The inflow is set by assigning the diffuser entrance conditions to cell 1=1. The outflow is set 
by extrapolating (with zero gradient) the flow from the inner boundary cell (I=L-1) to the 
boundary cell I = L. 
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(d) Conservation Laws (CYCEL'L) : 

Here all four conservation laws are integrated in time. There is an added equation of a "time- 
split" scheme for solving the conservation of H species equation. It is first solved without the 
recombination rate term (459), then the change in Z(l) due to recombination during that time 
step is added as DZZ (501 - 502). 

The energy equation includes the contribution of the formation energy indirectly. This equation 
has the same format as the adiabatic one (491), but when QDET.GT.O, the internal energy 
(505) is affected by changes in Z(l) and as a result the pressure P(l) (516) is also affected by 
progress in the recombination reaction. 



A. 2 Code Listing 



C$0PTI0NS LIST - 

IMPLICIT REAL^8(A-H,0 
C PROGRAM DETO 

C HFL -- HF/DF LASER DIFUSER 
COMMON B(52,26) 

1 .ENDB 
COMMON /AB/AC50) 
EQUIVALENCE (L,A(1)),(LL 

1 

2 

EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
COMMON 

1 






TRANSIENT FLOW. 



A(2) (T,A(3) ), (DT,A( A) ), (TMAX,A(5)), 
(TMUD,A(6)), (DTMUD,A(7)),(J0B,A(8)),(NERI,A(9) 
(JJJ,A(10)), (KEYM0N,A(11)),(NCYC,A(12)) 
(C0LELA,A(13)) 

(LAGEUL,A(14)) 

(UGAL,A(15)) 

(KEYEK,A(16)) 

(NCYCPR, A(17)) 

(STAB, A(18)), (DTBA,A(19)), (DTKOD,A(20) ) 
/M0NIT/NC1A(A),CASEAV(9),NF16(6), 

NMONUC 9 ) , NMONP ( ^ ) , NMONGC A ) 



(KDT,A(2D) 



NM0NR0(4),NM0NZ(9) 



HFLOOOl 

HFL0002 

HFL0003 

HFLOOOA 

HFL0005 

HFL0006 

HFL0007 

HFL0008 

HFL0009 

HFLOOlO 

HFLOOll 

HFL0012 

HFL0013 

HFL0014 

HFL0015 

HFL0016 

HFL0017 

HFL0018 






C 

C20 



DO 20 N=l,30 
NC1-^(N) = 0 
NMAT=26 

L=(LOCF(ENDB)-LOCF(B(l,l)))/NMAT 

L = 52 

LL=L-1 

NN=NMATXL 

DO 1 1=1, L 

DO 1 11=1, NMAT 

B(I,II)=0. 



HFL0020 

HFL0021 

HFL0022 

HFL0023 

HFL002A 

HFL0025 

HFL0026 

HFL0027 

HFL0028 

HFL0029 



CALL 


MAIN0(L,B(1 




1 




(] 


L, 


2 


),BC 


L, 


3) 


,B(] 


L, A 


),B(1 


. 5), 


HFL0030 


1 


B(l, 


6 


) 


,B( 


1, 


. 7 


) 


,B(1. 


• 8 


), 


B(l; 


• 9) 


,B(1, 


10), 


HFL0031 


2 


B(l, 


11 


) 


,B( 


1, 


. 12 


) 


,B(1; 


• 13 


), 


B(l; 


• lA) 


,B(1, 


15), 


HFL0032 


3 


B(l, 


16 


) 


,B( 


1, 


• 17 


) 


,B(1; 


• 18 


), 


B(l; 


• 19) 


,B(1, 


20), 


HFL0033 


A 


B(l, 


21 


) 


,B( 


1, 


• 22 


) 


,B(1; 


• 23 


), 


B(l; 


• 2A) 


,B(1, 


25), 


HFL003A 


5 


B(l, 


26 


) 


) 






















HFL0035 



STOP 



SUBROUTINE MAINO 

(L,X,U,P,RO,G,E,DU,DP, DRO, DG,DXSI,MIN, 
US,PS,UIDOT,PIDOT, 

FIMZ,ZMDOT, 

TENA, FIRO,FIM,FIE,GIP, VOL,Z,DZ) 

IMPLICIT REAL5(8(A-H,0-Z,$) 

DIMENSION X(L),U(L),P(L),RO(L),G(L),E(L),DU(L),DP(L),DRO(L), 
DG(L),DXSI(L),MIN(L), 

US(L),PS(L),UIDOT(L),PIDOT(L) 

,TENA(L), FIRO(L), FIM(L),FIE(L) 
,GIP(L),VOL(L),Z(L),DZ(L) 

,FIMZ(L),ZMDOT(L) 

COMMON /AB/A(50) 

EQUIVALENCE ( L L , A( 2 ) ) , ( T , A( 3 ) ) , ( DT , A( 9 ) ) , ( TMAX, A( 5 ) ) , 

(TMUD, A(6) ), (DTMUD, A(7)), ( J0B,A(8)), (NERI, A(9)), 

( JJJ, A(10) ), (KEYMON, A(1D), (NCYC,A(12) ) 
(LAGEUL,A(1A)) 

(NCYCPR,A(17)) 

(STAB,A(18)), (DTBA,A(19)), (DTKOD,A(20)),(KDT,A(2D) 



EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
COMMON /TOT/AMTOT, ETOT, EKTOT, EPTOT, TENTOT 



HFL0036 
HFL0037 
HFL d 038 
HFL0039 
HFLOOAO 
HFLOO-^l 
HFL00A2 
HFL00A3 
HFL0099 
HFL0095 
HFL00A6 
HFL0047 
HFLOO‘^8 
HFLOO‘^9 
HFL0050 
HFL0051 
HFL0052 
HFL0053 
HFL0054 
HFL0055 
HFL0056 
HFL0057 






T = 0 . 

NCYC=0 
JJJ = 0 

CALL NETUNM 
DELT=DT 
CALL BEGIN 



1 

2 

3 

( 

1 

2 



CALL SAFAE 



(L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN, 

US,PS,UIDOT,PIDOT, 

FIMZ,ZMDOT, 

TENA,FIRO,FIM, FIE, GIP, VOL , Z, DZ) 

(L,X,U,P,RO,G,E, DU, DP,DRO,DG,DXSI,MIN, 
US,PS,UIDOT,PIDOT, 

FIMZ,ZMDOT, 



HFL0059 

HFL0060 

HFL0061 

HFL0062 

HFL0063 

HFL0064 

HFL0065 

HFL0066 

HFL0067 

HFL0068 

HFL0069 

HFL0070 

HFL0071 

HFL0072 






LE: HFL 



SCRIPT A1 



3 TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ) 

1 NCYC=NCYC+1 

C TIME STEP CONTROL . 

DT = DTBA ■ 

IFCDT.GT.l . IDOXDTKOD.AND.DTKOD.NE. 0 . ) DT=1 . IDO^DTKOD 
IF(NCYC.EQ.Z) DT=DT/10. 

IF (NCYC.EQ.l) DT=0. 

IFCDT.EQ.O.) GO TO 11 
NHAD=((TMUD-T)/DT-1 »D-10) 

IFCNHAD.GE. 10) GO TO 11 
DT=(TMUD-T)/DFL0AT(NHAD+1) 

11 CONTINUE 
T=T+DT 

IF((NCYC/NCYCPR)xnCYCPR.NE.NCYC.AND.NCYC.GT.NCYCPR) go to 33 
PRINT 10, NCYC,T, DT,KDT 

10 FORMATCIX, *NCYC=M4,3X, ’ T= » , D1 1 » 4 , 3X, ’ DT= • , D1 1 . ^ , 3X, ’KDT=M4 
33 CONTINUE 

DTBA=DTMUD 
KDT = 0 
NERI=1 

IF (DABS(T-TMUD) .LT.l .D-8) NERI=0 
CALL CYCEUL 

1 (L,X,U,P,RO,G, E, DU, DP, DRO, DG, DXSI,MIN, 

2 US,PS,UIDOT,PIDOT, 

X FIMZ,ZMDOT, 

3 TENA,FIRO,FIM, FI E, GIP, VOL , Z, DZ) 

CALL SAFAE 

1 (L,X,U,P,RO,G, E,DU,DP,DRO,DG,DXSI,MIN, 

2 US,PS,UIDOT,PIDOT, 

FIMZ,ZMDOT, 

3 TENA, FIRO,FIM, FIE, GIP, VOL , Z, DZ) 

IF (NERI.NE.O) GO TO 2 

CALL PRINT 

1 (L,X,U,P,RO,G, E,DU,DP,DRO, DG,DXSI,MIN, 

2 US,PS,UIDOT,PIDOT, 

X FIMZ,ZMDOT, 

3 TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ) 

IF (DABS(T-TMUD) .LT.l.D-8) TMUD=TMUD+DTMUD 

2 CONTINUE 
DTKOD=DT 

IF (T.LT.TMAX-1 .D-8) GO TO 1 

RETURN 

END 

SUBROUTINE-NETURM 

IMPLICIT REAL^8(A-H,0-Z,$) 

COMMON /AB/AC50) 

EQUIVALENCE (L,A(D) 

EQUIVALENCE ( L L , A( 2 ) ) , ( T , A( 3 ) ) , ( DT, A( 4 ) ) , ( TMAX, A( 5 ) ) , 

1 (TMUD,A(6)), (DTMUD, A(7)), ( JOB, A(8) ), (NERI,A(9) ), 

2 ( JJJ,A(10)) , (KEYM0N,A(1D), (NCYC,A(12)) 

EQUIVALENCE (COL ELA, A( 13) ) 

EQUIVALENCE ( L AGEUL , A( 14 ) ) 

EQUIVALENCE ( KEYEK, A( 16 ) ) 

EQUIVALENCE ( NCYCPR, A( 17 ) ) 

EQUIVALENCE ( STAB , A ( 18 ) ) , ( DTBA , A ( 19 ) ) , ( DTKOD, A ( 20 ) ) , ( KDT , A ( 21 ) 
COMMON/DETO/QDET,TC, RATE, PCJDET, RCJDET, UCJDET, DCJDET, PODET, ROO 
C0MM0N/DIFFUS/U2,P2,R02, ARW 

COMMON /DRAW/GODELX, GODELY, UMIN, UMAX, PMIN, PMAX, ROMIN, ROMAX 
1 ,XMIN,XMAX,SMIN,SMAX, I VERSA 

COMMON /GAM/GAMA, NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10, Gil 
1 ,G12,G13,G14,G15,G16,G17,G18,G19,G20 

REALX8 NG 
REALMS MU2 

NAMELIST /IN/LIN, GAMA, DT,TMUD, DTMUD, TMAX, 

1 GODELX, GODELY, UMIN, UMAX, PMIN, PMAX, ROMIN, ROMAX, 

2 SMIN,SMAX, IVERSA,KEYMON, COL ELA, STAB 

3 ,LAGEUL, KEYEK 

4 ,QDET,TC,RATE 

LIN = L 

LAGEUL=2 

NCYCPR=20 



HFL0073 
HFL0074 
HFL0075 
HFL0076 
HFL0077 
HFL0078 
HFL0079 
HFL0080 
HFL0081 
HFL0082 
HFL0083 
HFL0084 
HFL0085 
HFL0086 
HFL0087 
HFL0088 
HFL0089 
HFL0090 
HFL0091 
HFL0092 
HFL0093 
HFL0094 
HFL0095 
HFL0096 
HFL0097 
HFL0098 
HFL0099 
HFLOlOO 
HFLOlOl 
HFL0102 
HFL0103 
HFL0104 
HFL0105 
HFL0106 
HFL0107 
HFL0108 
HFL0109 
HFLOllO 
HFLOlll 
HFL0112 
HFL0113 
HFL0114 
HFL0115 
"FFLUIT^- 
HFL0117 
HFL0118 
HFL0119 
HFL0120 
HFL0121 
HFL0122 
HFL0123 
HFL0124 
HFL0125 
HFL0126 
HFL0127 
HFL0128 
HFL0129 
HFL0130 
HFL0131 
HFL0132 
HFL0133 
HFL0134 
HFL0135 
HFL0136 
HFL0137 
HFL0138 
HFL0139 
HFL0140 
XXXXHFL0141 
HFL0142 
HFL0143 
HFL0144 



) 

DET 
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oooooo o oo 



HFL 



C 



SCRIPT A1 



KEYEK=1 




HFL0145 


TMUD=0. 




HFL0146 


DTMUD=0.2D0 




HFL0147 


TMAX=2.0D0 




HFL0148 


DT=2.5D-3 




HFL0149 


GAMA=1 .54D0 




HFL0150 


KEYM0N=1 




HFL0151 


STAB=0 .6D0 




HFL0152 


RATE=1 .5D0 




HFL0153 


G0DELX=16 . 




HFL0154 


GODELY=20. 




HFL0155 


IVERSA=100 




HFL0156 


UMIN=0. 




HFL0157 


UMAX= 3. DO 




HFL0158 


PMIN=0. 




HFL0159 


PMAX=1 .D-3 




HFL0160 


ROMIN=0 . 




HFL0161 


R0MAX=1 .D-3 




HFL0162 


SMIN=0. 




HFL0163 


SMAX=PMAX/ROMAX^^GAMA 




HFL0164 


COLELA=0. 




HFL0165 


READ IN 




HFL0166 

HFL0167 


PRINT IN 




HFL0168 

HFL0169 


GG=2.^GAMA/(GAMA-1 . ) 




HFL0170 


NG = GG 




HFL0171 


0 CONTINUE 




HFL0172 


MU2=(GAMA-1 . )/(GAMA+l . ) 




HFL0173 


G1=GAMA-1 . 




HFL0174 


G2=l .-MU2 




HFL0175 


G3 = 2./(3. ^GAMA-1 . ) 




HFL0176 


G4=(GAMA+1 . )/2. 




HFL0177 


G5=0 .5^(3.^GAMA-1 . )/(GAMA+l . ) 




HFL0178 


G6=(GAMA+1 . )/(2.^GAMA) 




HFL0179 


G7=2. /(GAMA-1 . ) 




HFL0180 


G8=( GAMA-1 . )/(2.^GAMA) 




HFL0181 


G9=(GAMA+1 . )/(4 .^GAMA) 




HFL0182 


G10=1./GAMA 




HFL0183 


G11=(GAMA+1. )/4. 




HFL0184 


G12=GAMA/(GAMA-1 . ) 




HFL0185 


G1 3=0. 5^( GAMA-3. )/(GAMA+l . ) 




HFL0186 


G14=0.5X(3.^GAMA-5. )/(GAMA+l . ) 




HFL0187 


G17=GAMA+1. 




HFL0188 


G0DELX=G0DELX/2.54 




HFL0189 


G0DELY=G0DELY/2.54 




HFL0190 


CALL NAMPLT(IVERSA) 




HFL0191 


CALL LIMITCIOOO.) 




HFL0192 


CALL PLOT(0.,0.5.-3) 




HFL0193 


DATA FOR DIFFUSER ENTRY CONDITIONS. TEST III/2, REF.(M-l). 


HFL0194 


UNITS ARE IN M,K,MILISEC. 




HFL0195 


FH2 IS THE DIFFUSER ENTRY MOLAR FRACTION OF 


H ATOMS. 


HFL0196 


FH2=0.091D0 




HFL0197 


U2=2.3D0 




HFL0198 


AMDOT=l .499D-5 




HFL0199 


A2=l .D0^7 .DO^O. 0254^^2 




HFL0200 


R02=AMD0T/(A2^U2) 




HFL0201 


TEMP0=1400.D0 




HFL0202 


W2=7 .27D0 




HFL0203 


ARW=8.31441D-3/W2 




HFL0204 


DELTAQ=435. 783/2. DO 




HFL0205 


QDET=FH2^DELTAQ/W2 




HFL0206 


AA=U2/DSQRT(GAMA^ARW^TEMP0) 




HFL0207 


RATE=1 .ADO^l .D7^FH2/W2^^2 




HFL0208 


RATE=RATE^1 .D2 




HFL0209 


EM2=AA/DSQRT(1 .-AA^^2/G7) 




HFL0210 


GOREM=(l .D0+EM2^^2/G7) 




HFL0211 


TRATI0=2 . D0^(GAMA+1 . DO )^EM2^^2^G0REM/ ( 1 . 


D0+GAMAXEM2^J«2)XS(2 


HFL0212 


HH0=GAMA^ARW^TEMP0/(GAMA-1 . DO) 




HFL0213 


HH1=HHO/TRATIO 




HFL0214 


CHOKE=QDET/(HH1-HHO) 




HFL0215 


TEMP2=TEMPO/GOREM 




HFL0216 
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P2 = ARWB^R02?^TEMP2 
RATIO=2.5D0/0.5D0 

EM0UT = EM2^RATI05()^( (GAMA-1 . DO)/ 2. DO) 

DO 250 1=1,100 
EM0UT1=EM0UT 

G0REM = RATI0B^EM0UT)^(1 .DO + 0 .5D0)^( GAMA-1 .DO) )^EM2B(X2) 
G0REM = G0REM)^^(2 . D0?^MU2) 

EMOUT = DSQRT( (2 . DO/ ( GAMA-1 . DO ) ) )( ( GOREM-1 . DO ) ) 
DM=DABS(EM0UT-EM0UT1) 

IFCDM.lt. 1 .D-9) GO TO 251 

250 CONTINUE 
CALL S0F(251) 

251 CONTINUE 
PRINT 201 

201 F0RMATC//1X, -DIFFUSER ENTRY DATA*) 

PRINT 202,U2,P2,R02,FH2 

202 FORMATC/IX, ‘ U2 , P2, R02 , FH2= ’ , AD16 . 5 ) 

PRINT 203, AMD0T,W2, TEMPO, TEMP2, EM2 

203 FORMAT (/IX, * AMDOT , W2, TEMPO , TEMP2, EM2= • , 5D16 . 5) 
PRINT 20A,RATE,QDET,HH0,HH1,CH0KE 

20 A FORMATC/IX, - RATE, QDET , HHO , HHl , CHOKE= ’ , 5D16 . 5/ ) 
PRINT 205, RATIO, EMOUT 

205 FORMATCIX, -EXIT AREA RATIO AND MACH NUMBER RATIO 
RETURN 
END 



subrouYine'be^in 

1 (L,X,U,P,RO,G, E,DU, DP,DRO,DG,DXSI,MIN, 

2 US,PS,UIDOT,PIDOT, 

FIMZ,ZMDOT, 

3 TENA,FIRO,FIM, FI E, GIP, VOL , Z, DZ) 

IMPLICIT REAU^8(A-H,0-Z,$) 

DIMENSION X(L),U(L),P(L),RO(L),G(L),E(L),DU(L),DP(L),DRO(L). 

1 DG(L),DXSI(L),MIN(L), 

2 US(L),PS(L),UIDOT(L),PIDOT(L) 

3 ,TENA(L),FIRO(L), FIM(L),FIE(L) 

A ,GIP(L), VOL(L),Z(L),DZ(L) 

5 ,FIMZ(L),ZMDOT(L) 

COMMON /AB/AC50) 

COMMON /GAM/GAMA, NG,MU2,G1,G2,G3,GA,G5,G6,G7,G8,G9,G10,G11 
1 ,G12,G13,G1A,G15,G16,G17,G18,G19,G20 

NG 



HFL0217 

HFL0218 

HFL0219 

HFL0220 

HFL0221 

^)^(0 .5D0/MU2)/EM2HFL0222 
HFL0223 
HFL022A 
HFL0225 
HFL0226 
HFL0227 
HFL0228 
HFL0229 
HFL0230 
HFL0231 
HFL0232 
HFL0233 
HFL023A 
HFL0235 
HFL0236 
HFL0237 
HFL0238 
HFL0239 
HFL02A0 
JflLQZfU, 



,EM0UT=-,2D16.5) 



REAU^8 
REAU^8 MU2 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 
EQUIVALENCE 



(LL,A(2)) 

(LAGEUL,A(1A)) 

(UGAL, A(15)) 

(STAB, A(18)), (DTBA,A(19)), (DTK0D,A(20)), (KDT, A(2D) 
COMMON/DETO/QDET,TC,RATE,PCJDET,RCJDET,UCJDET,DCJDET,PODET,ROODET 
COMMON /DRAW/GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX 
1 ,XMIN,XMAX,SMIN,SMAX, I VERSA 



HFL02A2 

HFL02A3 

HFL02AA 

HFL02A5 

HFL02A6 

HFL02A7 

HFL02A8 

HFL02A9 

HFL0250 

HFL0251 

HFL0252 

HFL0253 

HFL025A 

HFL0255 

HFL0256 

HFL0257 

HFL0258 

HFL0259 

HFL0260 

HFL0261 

HFL0262 

HFL0263 

HFL026A 

HFL0265 






DTBA=0. 


HFL0267 


DTK0D=0, 


HFL0268 


KDT = 0 


HFL0269 


P0=l.D-7 


HFL0270 


RHOO=l.D-7 


HFL0271 


Pl=0.1 


HFL0272 


P1=A . 


HFL0273 


RH01=0.125 


HFL027A 


RH01=A. 


HFL0275 


UGAL=0. 


HFL0276 


U0=1.08D0 


HFL0277 


U1 = 0. 


HFL0278 


UO=UO-UGAL 


HFL0279 


U1=U1-UGAL 


HFL0280 


UMIN=UMIN-UGAL 


HFL0281 


UMAX=UMAX-UGAL 


HFL0282 


X0=0.625D0 


HFL0283 


X1=3.125D0 


HFL028A 


XMIN=XO 


HFL0285 


XMAX=X1 


HFL0286 


DX=(Xl-X0)/(L-2. ) 


HFL0287 


DO 1 1=2, L 


HFL0288 
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X(I)=X0+(I“2. )^DX 


HFL0289 


1 


CONTINUE 


HFL0290 




X(L)=X1 


HFL0291 




XD=1 .DIO 


HFL0292 




DO 2 1=2, LL 


HFL0293 




U(I)=UO/X(I) 


HFL029A 




P(I)=PO 


HFL0295 




RO(I)=RHOO 


HFL0296 




Z(I)=0. 


HFL0297 




IF (X(I)-XD.LT.-l.D-8) GO TO 21 


HFL0298 




U(I)=U1 


HFL0299 




P(I)=P1 


HFL0300 




R0(I)=RH01 


HFL0301 




Z(I)=1.D0 


HFL0302 


21 


CONTINUE 


HFL0303 




GO TO (31,32), LAGEUL 


HFL030A 


31 


CONTINUE 


HFL0305 




E(I)=P(I)/((GAMA-1 . )^RO(I))+0 .5^U(I)^^2+Z(I)^QDET 


HFL0306 




GO TO 30 


HFL0307 


32 


CONTINUE 


HFL0308 




E(I)=P(I)/(GAMA-1 . )+0.5^RO(I)^U(I)^^2+Z(I)^RO(I)^QDET 


HFL0309 


30 


CONTINUE 


HFL0310 




G(I)=DSQRT(GAMA)^P(I)^RO(I)) 


HFL0311 


2 


CONTINUE 


HFL0312 




DO 3 1=2, LL 


HFL0313 




DXSI(I)=(X(I+1)’X(I))^R0(I) 


HFL031A 




TENA(I)=RO(I)^U(I) 


HFL0315 1 




VOL(I)=(X(I+1)-X(I))^(X(I+1)+X(I))/2.DO 


HFL0316 1 


3 


CONTINUE 


HFL0317 1 




RETURN 


HFL0318 i 




END 


HFL0319 




DOUBLE PRECISION FUNCTION RATIO(X) 


HFL0320 




IMPLICIT REAL^8(A-H,0"Z,$) 


HFL0321 




IFCX.LE.l .D-8)RATIO=0. 


HFL0322 




IF(X.LE.1.D-8)RETURN 


HFL0323 




RATIO=1.DO/X 


HFL032A 




RETURN 


HFL0325 




END 


HFL0326 



IMPLICIT REAL5(8(A-H,0-Z,$) HFL0328 

CROSS=X HFL0329 

RETURN HFL0330 

END HFL0331 

UUHkUUUNt CYChUL HHU332 ' 

1 (L,X,U,P,RO,G,E,DU,DP,DRO,DG,DXSI,MIN, HFL0333 

2 US,PS,UIDOT,PIDOT, HFL033A 

FIMZ,ZMDOT, HFL0335 

3 TENA,FIRO,FIM,FIE,GIP,VOL,Z,DZ) HFL0336 

IMPLICIT REALX8(A-H,0-Z,$) HFL0337 

DIMENSION X( L), U ( L), P( L), RO( L), G( L),E(L), DU (L),DP(L), DROLL), HFL 0338 

1 DG(L),DXSI(L),MIN(L), HFL0339 

2 US(L),PS(L),UIDOT(L),PIDOT(L) HFL03A0 

3 ,TENA(L),FIRO(L),FIM(L),FIE(L) HFL03A1 

A ,GIP(L),VOL(L),Z(L),DZ(L) HFL03A2 

5 ,FIMZ(L),ZMDOT(L) HFL03A3 

COMMON /AB/AL50) HFL03A4 

EQUIVALENCE ( LL , A( 2) ) , (T, A( 3) ) , ( DT, A( A ) ) , (COL ELA, A( 13) ) HFL03A5 

EQUIVALENCE ( KEYEK, A( 16 ) ) HFL03A6 

EQUIVALENCE ( STAB, A( 18 ) ) , ( DTBA , A( 19 ) ) , ( DTKOD, A( 20 ) ) , ( KDT , A( 21 ) ) HFL03A7 
COMMON /GAM/GAMA, NG,MU2,G1,G2,G3,GA,G5,G6,G7,G8,G9,G10, Gil HFL03A8 

1 ,G12,G13,G1A,G15,G16,G17,G18,G19,G20 HFL03A9 

REAL5«8 NG HFL0350 

REAL5<8 MU2 HFL0351 

COMMON /TOT/AMTOT,ETOT,EKTOT,EPTOT,TENTOT HFL0352 

COMMON /AZOV/ISAFA,NORIMN,USAF,PSAF,ROSAF,GSAF,ESAF,DPSAF HFL 0353 

1 ,DXSIL,DXSIR HFL035A 

LOGICAL NORIMN HFL0355 

COMMON /STEPO/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR, HFL0356 

1 RSTARL,RSTARR,GSTARL,GSTARR,WL,WR, HFL0357 

2 CL,CR,CSTARL,CSTARR,UW(6 ) HFL0358 

3 ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR HFL0359 

LOGICAL HELEML,HELEMR HFL0360 
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COMMON /STEP1/DUIDT,DPIDT,DGIDTL,DGIDTR,DRIDTL,DRIDTR HFL0361 

2 ,ASTARL,ASTARR HFL0362 

3 ,RAT,SH HFL0363 

COMMON ;'GRAD5/DUDXIL,DPDXIL,DGDXIL,DRDXIL, HFL0364 

1 DUDXIR,DPDXIR,DGDXIR,DRDXIR HFL0365 

COMMON /FI/FIH1,FIH2,FIH3,UXN,PXN,GXN,R0XN HFL0366 

1 ,GIH HFL0367 

2 ,FIHA,ZMDOTL,ZMDOTR HFL0368 

COMMON/DETO/QDET,TC,RATE,PCJDET,RCJDET,UCJDET,DCJDET,PODET,ROODET HFL 0369 

C DATA K0TZ/7777777777B/ HFL0370 



TN = T 


HFL0372 


DT2=DT/2. 


HFL0373 


DO 3 1=1, L 


HFL037A 


MINO=MIN(I) . AND.KOTZ 


HFL0375 


MINCI)=SHIFT(MINO,30) 


HFL0376 


CONTINUE 


HFL0377 


UXN=0 . 


HFL0378 


PXN=0 . 


HFL0379 


GXN=0 . 


HFL0380 


ZN = 0 . 


HFL0381 


ROXN=0 . 


HFL0382 


DO 1 1=2, L 


HFL0383 


IM=I“1 


HFL038^ 


UXNM=UXN 


HFL0385 


PXNM=PXN 


HFL0386 


GXNM=GXN 


HFL0387 


ROXNM=ROXN 


HFL0388 


ZNM=ZN 


HFL0389 


UL=U(IM)+0.5^DU(IM) 


HFL0390 


PL = P(IM)+0.55^DP(IM) 


HFL0391 


ROL=RO(IM)+0.5^DRO(IM) 


HFL0392 


GL=G(IM) + 0 .5)^DG(IM) 


HFL0393 


CL=GL/ROL 


HFL0394 


ZL = Z(IM) + 0.5)(DZ(IM) 


HFL0395 


UR=U(I)-0.5XDU(I) 


HFL0396 


PR=P(I)-0.5XDP(I) 


HFL0397 


GR=G(I)-0.5^DG(I) 


HFL0398 


ROR=RO(I)-0.5XDRO(I) 


HFL0399 


CR=GR/ROR 


HFLOAOO 


ZR=Z(I)-0.5XDZ(I) 


HFLOAOl 


CALL RI£MAN(L,I,MIN) 


HFL0A02 


DUDXIL=DU(IM)/DXSI(IM) 


HFL0A03 


DPDXIL=DP(IM)/DXSI(IM) 


HFLO^OA 


DGDXIL=DG(IM)/DXSI(IM) 


HFL0A05 


DRDXIL=DRO(IM)/DXSI(IM) 


HFL0^06 


DUDXIR=DU(I)/DXSI(I) 


HFL0A07 


DPDXIR=DP(I)/DXSI(I) 


HFL0A08 


DGDXIR=DG(I)/DXSI(I) 


HFL0A09 


DRDXIR=DRO(I)/DXSI(I) 


HFLOAIO 


SH=CROSS(X(I)) 


HFLOAll 


RAT=RATIO(X(D) 


HFL0A12 


CALL MAGA(L,I,MIN) 


HFL0A13 


US(I)=USTAR 


HFLOAIA 


PS(I)=PSTAR 


HFL0A15 


UIDOT(I)=DUIDT 


HFL0A16 


PIDOT(I)=DPIDT 


HFL0417 


CALL FLUXE1(L,I,MIN) 


HFL0A18 


CALL FLUXE2CL, I,MIN) 


HFL0A19 


FIR0(I)=FIH1 


HFL0^20 


FIM (I)=FIH2 


HFL0A21 


FIE (I)=FIH3 


HFL0A22 


GIP(I)=GIH 


HFL0A23 


FIMZ(I)=FIHA 


HFL042^4 


DU(IM)=UXN-UXNM 


HFL0A25 


DP(IM)=PXN-PXNM 


HFL0A26 


DG(IM)=GXN-GXNM 


HFL0A27 


DRO(IM)=ROXN-ROXNM 


HFL0428 


CONTINUE 


HFL0A29 


AMT0T=0 . 


HFL0A30 


ETOT=0 . 


HFL0A31 


EKT0T=0 . 


HFL0A32 
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EPT0T=0. 

TENTOT=0. 

FI1=FIR0(2) 

FI2=FIM (2) * - 
FIA=FIMZ(2) 

GI2=GIP(2) 

SH=CR0SS(X(2)) 

DO 2 1=2, LL 
IP=I+1 
FIM1=FI1 
FIM2=FI2 
FIMA=FIA 
GIM2=GI2 
SHM=SH 
GI2=GIP(IP) 

SH=CROSS(X(IP)) 

DVOL=VOL(I) 

FIl=FIRO(IP) 

FI2=FIM (IP) 

FIA=FIMZ(IP) 

ZKODM = Z(I))^RO(I) 

DX=X(IP)-X(I) 

RO(I)=RO( I)-DT/DV0U^(SH^FI1-SHM^FIM1) 

TENA(I)=TENA(I)“DT/DV0LX(SH)^FI2-SHMXFIM2)-DT/DXX(GI2-GIM2) 

U(I)=TENA(I)/RO(I) 

IFCQDET.EQ.O. ) GO TO 2 

Z(I) = (ZKODM”(DT/DVOL))^(SH)^FIA-SHMxFIMA) )/RO(I) 

DZZl=Z(I)-l.DO 

IF(DZZ1.LT.1.D"6.AND.Z(I) .GE.O.) GO TO 7350 
IFCDZZl.LT.O.) GO TO 7352 
Z(I)=I .DO 

IFCDZZl .LT.O .2D0) GO TO 7350 

7352 CONTINUE 

IF(Z(I) .GT.O. ) GO TO 7353 
Z(I)=0. 

IF(DABS(Z(D) .LT.0.2D0) GO TO 7350 

7353 CONTINUE 

PRINT 7351, I,Z(I),P(I),RO(I), E(I),U(I),ZMDOT(I) 

7351 FORMATC/IX, 'CYCEUL . I , Z, P, RO, E, U, ZMDOT= M5, 6 D15 . 5/) 

CALL SOFC7350) 

7350 CONTINUE 

IF(Z(I) .GT. .IDl) Z(I)=.1D1 
2 CONTINUE 

IFCCOLELA.EQ.O. ) GO TO 201 
CALL DCOLE(L,X,U , DU ,MIN,1) 

CALL DCOLE(L,X,RO,DRO,MIN, A) 

201 CONTINUE 

CALL BD0K1(L,X,U , DU ,MIN,1) 

CALL BD0K1(L,X,R0,DR0,MIN, A) 

FI3=FIE (2) 

SH=CR0SS(X(2)) 

DO 6 1=2, LL 
IP=I+1 
SHM=SH 

SH=CROSS(X(IP)) 

DVOL=VOL(I) 

FIM3=FI3 
FI3=FIE (IP) 

DX=X(IP)-X(I) 

E(I)=E(I)-DT/DV0U^(SHXFI3-SHM^FIM3) 

C 

UAV=U(I) 

C 

ROAV=RO(I) 

C 

EP = E( I )-0 . 5?^R0AV?^UAV?(X2 
IF (QDET.EQ.O.) GO TO 6A 
TI=P(I) 

DZZ = -RATEXDT5^(RO(I)^Z(I))X)^2 

Z(I)=Z(I)+DZZ 

IF(Z(I) .LT.O. ) Z(I)=0. 

IF (Z(I) .LT.O.OI) Z(I)=0. 



HFL0A33 

HFL0A3A 

HFL0A35 

HFL0A36 

HFL0A37 

HFL0A38 

HFL0A39 

HFLOAAO 

HFLOAAl 

HFL0AA2 

HFL0AA3 

HFLOAAA 

HFL0AA5 

HFL0AA6 

HFL0AA7 

HFL0AA8 

HFL0AA9 

HFL0A50 

HFL0A51 

HFL0A52 

HFL0A53 

HFL0A5A 

HFL0A55 

HFL0A56 

HFL0A57 

HFL0A58 

HFL0A59 

HFL0A60 

HFL0A61 

HFL0A62 

HFL0A63 

HFL0A6A 

HFL0A65 

HFL0A66 

HFL0A67 

HFL0A68 

HFL0A69 

HFL0A70 

HFL0A71 

HFL0A72 

HFL0A73 

HFL0A7A 

HFL0A75 

HFL0A76 

HFL0A77 

HFL0A78 

HFL0A79 

HFL0A80 

HFL0A81 

HFL0A82 

HFL0A83 

HFL0A8A 

HFL0A85 

HFL0A86 

HFL0A87 

HFL0A88 

HFL0A89 

HFL0A90 

HFL0A91 

HFL0A92 

HFL0A93 

HFL0A9A 

HFL0A95 

HFL0A96 

HFL0A97 

HFL0A98 

HFL0A99 

HFL0500 

HFL0501 

HFL0502 

HFL0503 

HFL050A 
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EP=EP-Z(I)xRO(I)xQDET 


HFL0505 


6A 


CONTINUE 


HFL0506 




GO TO .(61,62), KEYEK 


HFL0507 


61 


CONTINUE 


HFL0508 




GO TO 60 


HFL0509 


62 


CONTINUE 


HFL0510 




DEK=(ROAVxDU( I)x^2+2 . XDROC I )XDU( I )XUAV)/2A . 


HFL0511 




EP=EP-DEK 


HFL0512 




GO TO 60 


HFL0513 


60 


CONTINUE 


HFL051A 




IFCEP.LE.O.) GO TO 7001 


HFL0515 




P(I)=G1XEP 


HFL0516 




G(I)=D5QRT(GAMAXP(I)XR0(I) ) 


HFL0517 




DM=RO(I)XDX 


HFL0518 




DXSI(I)=DM 


HFL0519 




DM=RO(I)xDVOL 


HFL0520 




ETOT=ETOT+E(I)XDVOL 


HFL0521 




EPTOT=EPTOT+EPXDVOL 


HFL0522 




AMTOT=AMTOT+DM 


HFL0523 




ETOT=ETOT+E(I)xDX 


HFL052A 




EPT0T=EPT0T+EPXDX 


HFL0525 




TENTOT=TENTOT+DXXTENA(I) 


HFL0526 




UPC=DABS(U(I) )+G(I)/RO(I) 


HFL0527 




DTI=STABxDX/UPC 


HFL0528 




IFCDTI .GT.DTBA) GO TO 69 


HFL0529 




DTBA=DTI 


HFL0530 




KDT = I 


HFL0531 


69 


CONTINUE 


HFL0532 


6 


CONTINUE 


HFL0533 




CALL DC0LE(L,X,Z,DZ,MIN,5) 


HFL053A 




IFCCOLELA.EQ.O. ) GO TO 601 


HFL0535 




CALL DC0LE(L,X,P,DP,MIN,2) 


HFL0536 




CALL DC0LE(L,X,G,DG,MIN,3) 


HFL0537 


601 


CONTINUE 


HFL0538 




CALL BD0K1(L,X,P,DP,MIN,2) 


HFL0539 




CALL BD0K1(L,X,G,DG,MIN,3) 


HFL05A0 




CALL BD0K1(L,X,Z,DZ,MIN,5) 


HFL05A1 



EKTOT=ETOT-EPTOT 

RETURN 

7001 CONTINUE 

PRINT 7101, I,ROAV,UAV;DRO(I),DU(I),E(I),DEK,EP,KEYEK 
7101 F0RMATC//1X, 'FROM CYCEUL . NEGATIVE EP . IN CELL I = M6// 

1 IX, * ROAV, UAV, DRO( I ) , DU( I) = * , AD18 . 8// 

2 IX, * E(I),DEK,EP,KEYEK=' ,3D18 .8,110//) 

CALL SOF(7001) 

RETURN 
FNn 



HFL05A2 
HFL05A3 
HFL05AA 
HFL05A5 
HFL05A6 
HFL05A7 
HFL05A8 
HFL05A9 
HFL0550 
■ HFLQ^5 , 1 . 



.DU(L),DP(L),DRO(L). 



SUBROUTINE SAFAE 

1 (L,X,U,P,RO,G, E, DU,DP,DRO,DG, DXSI,MIN, 

2 US,PS,UIDOT,PIDOT, 

^ FIMZ,ZMDOT, 

3 TENA,FIRO, FIM,FIE,GIP,VOL,Z,DZ) 

IMPLICIT REAU(8(A-H,0-Z,$) 

DIMENSION X(L),U(L),P(L),RO(L),G(L),E(L). 

1 DG(L),DXSI(L),MIN(L), 

2 US(L),PS(L),UIDOT(L),PIDOT(L) 

3 ,TENA(L),FIRO(L),FIM(L),FIE(L) 

A ,GIP(L), VOL(L),Z(L),DZCL) 

5 ,FIMZ(L),ZMDOT(L) 

COMMON /AB/AC50) 

EQUIVALENCE ( L L , A( 2 ) ) , ( T , A( 3 ) ) , ( DT , A( A ) ) , ( NCYC, A(12 ) ) 

COMMON /GAM/GAMA, NG, MU2, Gl, G2, G3, GA, G5, G6 ,G7,G8,G9,G10,G11 
1 ,G12,G13,G1A,G15,G16,G17,G18,G19,G20 

REALMS NG 
REALMS MU2 

COMMON/DETO/QDET,TC,RATE,PCJDET,RCJDET,UCJDET, DCJDET,PODET,ROODET 
COMMON/ DIFFUS/U2,P2,R02, ARM 

TN=T HFL0573 

C HFL057A 

C INFLON B.C. AT 1=2 HFL0575 

C HFL0576 



HFL0552 

HFL0553 

HFL055A 

HFL0555 

HFL0556 

HFL0557 

HFL0558 

HFL0559 

HFL0560 

HFL0561 

HFL0562 

HFL0563 

HFL056A 

HFL0565 

HFL0566 

HFL0567 

HFL0568 

HFL0569 

HFL0570 

HFL0571 
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U(1)=U2 


HFL0577 




P(1)=P2 


HFL0578 




R0(1)=R02 


HFL0579 




G(l)=DSQRT(GAMAxPa)?^ROa) ) 


HFL0580 




Z(l)=l .DO 


HFL0581 




DU(1)=0. 


HFL0582 




DP(1)=0. 


HFL0583 




DG(1)=0. 


HFL0584 




DRO(1)=0. 


HFL0585 




DXSI(1)=DXSI(2) 


HFL0586 


C 




HFL0587 


C OUTFLOW B.C. AT I=L 


HFL0588 


C 




HFL0589 




U(L)=U(LU 


HFL0590 




P(L)=P(LL) 


HFL0591 




RO(L)=RO(LL) 


HFL0592 




GCL)=DSQRT(GAMAB(P(L)B^RO(L)) 


HFL0593 




DU(L)=0. 


HFL0594 




DP(L)=0. 


HFL0595 




DG(L)=0. 


HFL0596 




DRO(L)=0. 


HFL0597 




DXSI(L)=DXSI(LL) 


HFL0598 


C 




HFL0599 




RETURN 


HFL0600 




END 


HFi n^m 




SUBROUTINE BDOKl ( L , X, V, DV , MIN, NV) 


HFL0602 




IMPLICIT REALB^3(A-H,0-Z,$) 


HFL0603 




DIMENSION X(L),V(L),DV(L),MIN(L) 


HFL0604 




COMMON /AB/AC50) 


HFL0605 




EQUIVALENCE ( L L , A( 2 ) ) , ( KEYMON, A( 1 1 ) ) 


HFL0606 




COMMON /DRAW/GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX 


HFL0607 




1 ,XMIN,XMAX,SMIN,SMAX,IVERSA 


HFL0608 




COMMON /M0NIT/NC1A(4),CASEAV(A),NF16C6), 


HFL0609 




1 NM0NU(4),NM0NP(^),NM0NG(A),NM0NR0(^) .NMONZ(A) 


HFL0610 




DIMENSION NM0NV(A,5) 


HFL0611 




EQUIVALENCE ( NMONVC 1 , 1) , NMONUC 1 ) ) 


HFL0612 


C 


DIMENSION NAMEVC5) 


HFL0613 


C 


DATA MAMEV/1HU,1HP,1HG,2HR0,1HZ/ 


HFL0614 




DATA EPS/1. D-15/ 


HFL0615 




c 


NV = 0 


HFL0617 


c 


DO 10 N=l,5 


HFL0618 


c 


IF (NAME.EQ.NAMEV(N))NV=N 


HFL0619 


CIO 


CONTINUE 


HFL0620 




GO TO (1,2, 3, 4, 5), NV 


HFL0621 


1 


AMIDA = (UMAX-UMIN)BCX2 


HFL0622 




GO TO 9 


HFL0623 


2 


AMIDA = (PMAX-PMIN)B(X2 


HFL0624 




GO TO 9 


HFL0625 


3 


AMI DA = ((UMAX-UMIN)x(ROMAX-ROMIN)) 5^X2 


HFL0626 




GO TO 9 


HFL0627 


A 


AMIDA = (R0MAX-R0MIN)XB^2 


HFL0628 




GO TO 9 


HFL0629 


5 


AMIDA=1 .DO 


HFL0630 




GO TO 9 


HFL0631 


9 


CONTINUE 


HFL0632 




AMIDA=AMIDA^EPSxx2 


HFL0633 




AMIDA=EPS 


HFL0634 




DO 29 1=2, LL 


HFL0635 




ICAT=0 


HFL0636 




VLEFT = V(I)-0.5D0B(DV(I) 


HFL0637 




VRIGHT=V(I)+0.5D0^DV(I) 


HFL0638 




VM=V(I-1) 


HFL0639 




VP=V(I+1) 


HFL0640 




SIGN=(VP-V(I) )X(V(I)-VM) 


HFL0641 




IF(SIGN.GT.-AMIDA) GO TO 22 


HFL0642 


21 


DV(I)=0. 


HFL0643 




ICAT=1 


HFL0644 




GO TO 20 


HFL0645 


22 


CONTINUE 


HFL0646 




SIGN = (VP-VM)B^DV(I) 


HFL0647 




IF(SIGN.GT.-AMIDA) GO TO 24 


HFL0648 



30 



ILE: 



HFL 





SCRIPT A1 




23 


DV(I)=0.5D0*(VP-VM) 


HFL0649 




VLEFT=V(I)-0.5D0XDV(I) 


HFL0650 




VRIGHT=V(I)+0 .5D0XDVCI) 


HFL0651 




ICAT=2 


HFL0652 


24 


SIGN=(VLEFT-VM)XDV(I) 


HFL0653 




IF(SIGN.GT.-AMIDA) GO TO 26 


HFL0654 


25 


VLEFT=VM 


HFL0655 




VRIGHT = 2. D0*V( D-VLEFT 


HFL0656 




DV(I)=VRIGHT-VLEFT 


HFL0657 




ICAT=3 


HFL0658 


26 


SIGN=(VP-VRIGHT)*DV( I) 


HFL0659 




IF (SIGN.GT. -AMIDA) GO TO 28 


HFL0660 


27 


VRIGHT=VP 


HFL0661 




VLEFT = 2. DOXVCD-VRIGHT 


HFL0662 




DV(I)=VRIGHT-VLEFT 


HFL0663 




ICAT=3 


HFL0664 


28 


IF(DABS(DV(I)).LE.0.5D0XDABS(VP-VM)) GO TO 31 


HFL0665 


30 


DV(I)=0.5D0X(VP-VM) 


HFL0666 




ICAT=4 


HFL0667 


31 


CONTINUE 


HFL0668 


20 


CONTINUE 


HFL0669 




IF (DV(I)xx2.GT. AMIDA) GO TO 40 


HFL0670 




DV(I)=0. 


HFL0671 


40 


CONTINUE 


HFL0672 


C 


IF (ICAT.GT.O) NMONV(ICAT,NV)=NMONV(ICAT,NV)+l 


HFL0673 


C 


IBYTE0=5 


HFL0674 


C 


MIN(I)=MIN(I) .OR.SHIFTCICAT, (NV+IBYTE0-2)X3) 


HFL0675 


29 


CONTINUE 


HFL0676 




RETURN 


HFL0677 




END _ . ___ . . . 


HFLQ67_8. 




SUBROUTINE DCOL E( L , X, V , DV , MIN , NV ) 


HFL0679 




IMPLICIT REALX8(A-H,0-Z,$) 


HFL0680 




DIMENSION X(L),V(L),DV(L),MIN(L) 


HFL0681 




COMMON /AB/AC50) 


HFL0682 




EQUIVALENCE (LL,A(2)) 


HFL0683 






DO 1 1=2, LL 

IM=I-1 

IP=I+1 

DV(I)=0 .5X(V(IP)-V(IM)) 
CONTINUE 
RETURN 
JEMU, 



HFL0685 

HFL0686 

HFL0687 

HFL0688 

HFL0689 

HFL0690 

HFL0691 



SUBROUTINE PRINT 

1 (L,X,U,P,RO,G, E,DU,DP, DRO, DG, DXSI , MIN, 

2 US,PS,UIDOT,PIDOT, 

X FIMZ,ZMDOT, 

3 TENA, FIRO, FIM,FIE,GIP,VOL,Z,DZ) 
IMPLICIT REALX8CA-H,0-Z,$) 

DIMENSION X(L),U(L),P(L),RO(L),G(L),E(L),DU(L) ,DP(L), 



DRO(L), 



1 DG(L),DXSI(L),MIN(L), 

2 US(L),PS(L),UIDOT(L),PIDOTCL) 

3 ,TENA(L), FIRO(L),FIM(L),FIE(L) 

4 ,GIP(L),VOL(L),Z(L),DZ(L) 

5 ,FIMZ(L),ZMDOT(L) 

COMMON /TOT/AMTOT, ETOT, EKTOT, EPTOT, TENTOT 

COMMON /STEPO/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR, 

1 RSTARL,RSTARR,GSTARL,GSTARR,WL,WR, 

2 CL,CR,CSTARL,CSTARR,UW(6) 

3 ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR 
LOGICAL HELEML,HELEMR 

COMMON /AB/AC50) 

EQUIVALENCE ( L L , A( 2) ) , ( T, A ( 3 ) ) , ( NCYC, A( 12) ) , ( DT , A ( 4) ) 

EQUIVALENCE (UGAL,A(15)) 

COMMON/DETO/QDET,TC,RATE,PCJDET,RCJDET,UCJDET,DCJDET,PODET,ROODET 

C0MM0N/DIFFUS/U2,P2,R02,ARW 

COMMON /GAM/GAMA, NG,MU2,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10, Gil 
1 ,G12,G13,G14,G15,G16,G17,G18,G19,G20 

REALX8 NG 
REALX8 MU2 

COMMON /M0NIT/NC14(4),CASEAV(4),NF16(6), 

1 NMONU ( 4 ) , NMONP ( 4 ) , NMONGC 4 ) , NMONRO ( 4 ) , NMONZC 4 ) 



HFL()692 

HFL0693 

HFL0694 

HFL0695 

HFL0696 

HFL0697 

HFL0698 

HFL0699 

HFL0700 

HFL0701 

HFL0702 

HFL0703 

HFL0704 

HFL0705 

HFL0706 

HFL0707 

HFL0708 

HFL0709 

HFL0710 

HFL0711 

HFL0712 

HFL0713 

HFL0714 

HFL0715 

HFL0716 

HFL0717 

HFL0718 

HFL0719 

HFL0720 



31 



E: HFL 



SCRIPT A1 



DIMENSION CASAV1C4) HFL0721 
LOGICAL FULLER HFL0722 
COMMON /PRNT/ROX,PX,UX,ZX HFL0723 



44 



11 



12 



222 



FULLPR=.TRUE. 

PRINT 1 
FORMATCIHI) 

PRINT 2, T,DT,NCYC 

FORMATClXaOX, 'RESULTS AT T= S D1 1 . 5 , 5X. ' DT= S D1 1 . 5 , 5X, »NCYC=’ 
1 15//) 

PRINT 3, AMTOT,ETOT,EKTOT,EPTOT,TENTOT 



’AMTOT=' 
'TENTOT= 
' I’ 



, D20 . 14.2X, 'ETOT, EKTOT, EPTOT=V 
SD21 .14//) 



3D22.14/ 



13 

131 

10 



FORMATCIX, 

1 IX, 

FORMATdX,' I’,’ X U P *, 

1 * RO G Z *, 

2 * DU DP DRO *, 

3 * DG * , * DZ * ) 

FORMATCIX,' *,* US PS 

1 * ZMDOT S' FIMZ ',' AMDOT ', 

2 ' AMDOTN ',' TEMP ',' ENTALP ', 

3 ' AMACH ',' ENTRO ') 

FORMAT(IX) 

IF (UGAL.NE.O.) PRINT 6, UGAL 

FORMATC/llX, 'INITIAL VELOCITY CORRESPONDS TO UGAL = ' , D15 . 6/ ) 

DO 10 1=1, L 

IF (MOD(I,10) .NE.l) GO TO 11 

PRINT 5 

PRINT 4 

PRINT 44 

PRINT 5 

CONTINUE 

PRINT 12, 1,X(I) ,U(I),P(I),RO(I) ,G(I),Z(I),DU(I),DP(I) ,DRO(I), 

1 DGa),D2(I) 

F0RMAT(1X,I3,6D12.5,5D11 .4) 

ENTRO=P(I)/RO(I)xxGAMA 
IF( .NOT. FULLER) GO TO 131 
IF(I.EQ.l) GO TO 131 
IM=I-1 

UL = U(IM) + 0.5^eDU(IM) 

PL=P(IM)+0 .S^DPCIM) 

ROL=ROCIM)+0.5XDRO(IM) 

GL=G(IM) + 0.5^(DG(IM) 

CL=GL/ROL 

ZL=Z(IM)+0.5^DZ(IM) 

UR=U(I)-0.5XDUCI) 

PR=P(I)-0.5^DP(I) 

GR=G(I)-0.5^DG(I) 

ROR=RO(I)-0.5XDRO(I) 

CR=GR/ROR 

ZR=Z(I)-0.5^DZ(I) 

CALL RIEMANCL,I,MIN) 

CALL FLUXE1CL,I,MIN) 

XI=X(I) 

IF(I.EQ.L) GO TO 222 
UX=U(I) 

PX=P(I) 

ROX=RO(I) 

ZX=Z(I) 

IP=I+1 

XI=0.5D0X(X(I)+X(IP)) 

CONTINUE 

AMACH = UX/DSQRT(GAMA^cpX/ROX) 

AMDOT=ROXxUXxCROSS(XI) 

IFCI . EQ.2)AMDOTO=AMDOT 
AMDOTN=AMDOT/AMDOTO 

ENTAL P = ( GAMA/ C GAMA-1 .DO) )XPX/ROX+0 .5D0XUX^3(2 
TEMP = PX/(ROX5(ARW) 

PRINT 13,US(I),PS(I), 

1 ZMDOT (I), FIMZ (I), AMDOT, AMDOTN, TEMP, ENTALP, AMACH, ENTRO 



F0RMATC4X, 12X,5D12.5,6D11 

CONTINUE 

CONTINUE 



4) 



HFL0725 

HFL0726 

HFL0727 

HFL0728 

HFL0729 

HFL0730 

HFL0731 

HFL0732 

HFL0733 

HFL0734 

HFL0735 

HFL0736 

HFL0737 

HFL0738 

HFL0739 

HFL0740 

HFL0741 

HFL0742 

HFL0743 

HFL0744 

HFL0745 

HFL0746 

HFL0747 

HFL0748 

HFL0749 

HFL0750 

HFL0751 

HFL0752 

HFL0753 

HFL0754 

HFL0755 

HFL0756 

HFL0757 

HFL0758 

HFL0759 

HFL0760 

HFL0761 

HFL0762 

HFL0763 

HFL0764 

HFL0765 

HFL0766 

HFL0767 

HFL0768 

HFL0769 

HFL0770 

HFL0771 

HFL0772 

HFL0773 

HFL0774 

HFL0775 

HFL0776 

HFL0777 

HFL0778 

HFL0779 

HFL0780 

HFL0781 

HFL0782 

HFL0783 

HFL0784 

HFL0785 

HFL0786 

HFL0787 

HFL0788 

HFL0789 

HFL0790 

HFL0791 

HFL0792 
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C 
C 

c 

c 

CAO 

c 

C30 

C 

C31 

C 

c 

C301 

C 

C 

C32 

C 

C 

C 

C33 

C 

C 

C 

C 

C 



JOB 



,NE. 0) CASAV1(I)=CASEAV(I)/NC1A(I) 



STATISTICS S3X, 10 

IN RIEMAN SOLVER NC1A(NCASE)= 



1 



STATISTICS 
DO AO 1=1, A 
CASAV1(I)=0 
IF (NCl-A(I) 

CONTINUE 
PRINT 30 

FORMAT (/// IX, 10(1H^^) ,3X, 'JOB 
PRINT 31, (NC1A(I),I=1, A) 

FORMATCIX, 'NO. OF VARIOUS CASES 
AIIO) 

PRINT 301, (CASAVl(I), 1=1, A) 

FORMATC/IX, 'AVERAGE NUMBER OF ITERATIONS IN RIEMAN SOLVER', 
IX, ' CASAV1(NCASE)=',A(F6.2,AX)) 

PRINT 32, (NF16(I), 1=1,6) 

FORMATC/IX, 'NO. OF VARIOUS FLUX CASES NF16 ( NFLUX) = ' , 6 1 10 ) 

ICAT0=A 

PRINT 33, (NMONU(I) , I=1,ICAT0), (NMONP(I), 1=1, ICATO), 
(NMONG(I), 1=1, ICATO), (NMONRO(I), 1=1, ICATO) 
F0RMAT(/1X, 'NO. OF MONOTONICITY INTERVENTIONS FOR EACH VAR.' 
IX, 'IN EACH CATEGORY. '/ 

IX, 'NMONU dCAT) = ' , AIIO/ 

IX, 'NMONP (ICAT)=',AI10/ 

IX, 'NMONG (ICAT)=',AI10/ 

IX, 'NMONRO(ICAT)=', AIIO/) 

RETURN 



HFL0793 

HFL079A 

HFL0795 

HFL0796 

HFL0797 

HFL0793 

HFL0799 

HFL0800 

HFL0801 

HFL0302 

HFL0805 

HFL080A 

HFL0805 

HFL0806 

HFL0807 

HFL0808 

HFL0809 

HFL0810 

HFL0811 

HFL0812 

HFL0813 

HFL081A 

HFL0815 

HFL0816 

HFL0817 





SUBROUTINE SOF(ISTOP) 


HFL0819 




IMPLICIT REALX8(A-H,0-Z,$) 


HFL0820 




PRINT 1, ISTOP 


HFL0821 


1 


FORMATC lX,3H^^x^,3X, 'SOF' , 16 , 3X, 3H^^^^/ ) 


HFL0822 




XX=-1 . DO 


HFL0823 




YY=DSQRT(XX) 


HFL082A 


c 


CALL SYSTEMC51, 16H SUBROUTINE TREE) 


HFL0825 




STOP 


HFL0826 







HFL0827 



SUBROUTINE RI EMANC L , I , MI N) 

IMPLICIT REALX8(A-H,0-Z,$) 

DIMENSION MIN(L) 

COMMON /STEPO/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR, 

1 RSTARL,RSTARR,GSTARL,GSTARR,WL,WR, 

2 CL,CR,CSTARL,CSTARR,UW(6) 

3 ,ZL,ZR,ZSTARL,2STARR,NFLUX,HELEML,HELEMR 
LOGICAL HELEML,HELEMR 

COMMON /STEPl/DUIDT, DPIDT, DGIDTL, DGIDTR, DRIDTL, DRIDTR 

2 , ASTARL, ASTARR 

3 ,RAT,SH 

COMMON /DRAW/GODELX,GODELY,UMIN,UMAX,PMIN,PMAX,ROMIN,ROMAX 
,XMIN,XMAX,SMIN,SMAX,IVERSA 



COMMON /GAM/GAMA, NG,MU2,G1,G2,G3,GA,G5,G6,G7,G8,G9,G10,G11 



,G12,G13,G1A,G15,G16,G17,G18,G19,G20 

REALMS NG 
REAU^S MU2 
COMMON /AB/AC50) 

COMMON /M0NIT/NC1A(A),CASEAV(A) ,NF16(6), 

1 NMONUC A),NMONP(A),NMONG(A), NMONRO(A),NMONZ(A) 

DATA NMAX/63/ HFL08A9 

DATA EPS/1. D-8/ HFL0850 



HFL0828 

HFL0829 

HFL0830 

HFL0831 

HFL0832 

HFL0833 

HFL083A 

HFL0835 

HFL0836 

HFL0837 

HFL0838 

HFL0839 

HFL08A0 

HFL08A1 

HFL08A2 

HFL08A3 

HFL08AA 

HFL08A5 

HFL08A6 

HFL08A7 



7211 



IFCPL.GT.O. .AND.PR.GT.O.) GO TO 7201 


HFL0852 


PRINT 7211, I, PL, PR 


HFL0853 


FORMATC' FROM RIEMAN. I , PL , PR= ' , 16 , 2D16 . 6/ ) 


HFL085A 


CALL S0FC765A) 


HFL0855 


CONTINUE 


HFL0856 


UW(6)=1 .D20 


HFL0857 


WL = 0 . 


HFL0858 


NR = 0 . 


HFL0859 


ZETAL=PL^^G8 


HFL0860 


ZETAR = PR)^XG8 


H FL 0 8 6 1 


CLG=CL/GAMA 


HFL0862 


CRG=CR/GAMA 


HFL0863 


ZSTARL=ZL 


HFL0864 



33 
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ZSTARR=ZR 


HFL0865 




IF (ZETAL .LT.ZETAR) GO TO 102 


HFL0866 


C LEFT PRESSURE IS HIGHER 


HFL0867 


101 


CONTINUE 


HFL0868 




EVERR=(PL“PR)/PR 


HFL0869 




USR = UR+CRG3(EVERR/DSQRT(1 .+G6XEVERR) 


HFL0870 




SR=USR 


HFL0871 




UEL = UL-G7 5^CLX(ZETAR-ZETAL)/ZETAL 


HFL0872 




SL=UEL 


HFL0873 




NL = 2 


HFL0874 




NR = 2 


HFL0875 




IF (USR.GE,UL) NL=1 


HFL0876 




IF (UEL.LE.UR) NR=1 


HFL0877 




IF (DABS(EVERR) .LT.EPS) GO TO 100 


HFL0878 




IF (NL.EQ.2.AND.NR.EQ.1) GO TO 7001 


HFL0879 




GO TO 100 


HFL0880 


C RIGHT PRESSURE IS HIGHER 


HFL0881 


102 


CONTINUE 


HFL0882 




EVERL=(PR-PL)/PL 


HFL0883 




USL = UL-CLGXEVERL/DSQRT(1 .+G63(EVERL) 


HFL0884 




SL=USL 


HFL0885 




UER = UR+G7)^CRX( ZETAL -ZETAR)/ZETAR 


HFL0886 




SR=UER 


HFL0887 




NL=2 


HFL0888 




NR = 2 


HFL0889 




IF (UER.GE.UL) NL=1 


HFL0890 




IF (USL.LE.UR) NR=1 


HFL0891 




IF (DABS(EVERL) .LT.EPS) GO TO 100 


HFL0892 




IF (NL.EQ.1.AND.NR.EQ.2) GO TO 7001 


HFL0893 




GO TO 100 


HFL0894 


100 


CONTINUE 


HFL0895 




IF (NL.EQ.1.AND.NR.EQ.2) NCASE=1 


HFL0896 




IF (NL.EQ.2.AND.NR.EQ.2) NCASE=2 


HFL0897 




IF (NL.EQ.2.AND.NR.EQ.1) NCASE=3 


HFL0898 




IF (NL.EQ.l.AND.NR.EQ.l) NCASE=4 


HFL0899 




IF(DABS(PL-PR)+DABS(UL-UR) . LT . EPSX( PMAX-UMI N ) ) NCASE=4 


HFL0900 




UMIDA=EPSXDMAX1(CL,CR) 


HFL0901 




DUDZL=-G7XCL/ZETAL 


HFL0902 




DUDZR= G7XCR/ZETAR 


HFL0903 




ZETA=(-(UR-UL)+ZETAR^DUDZR-ZETALXDUDZL)/(DUDZR-DUDZL) 


HFL0904 




IF (ZETA.LE.O. ) GO TO 7002 


HFL0905 




N = 0 


HFL0906 




GO TO (1,2, 3,4), NCASE 


HFL0907 


C THE CASE ES 


HFL0908 


1 


ITYPE=NCASE 


HFL0909 




HELEML=. FALSE. 


HFL0910 




HELEMR= .TRUE. 


HFL0911 


11 


N = N+1 


HFL0912 




IF (N.GT.NMAX) GO TO 7003 


HFL0913 




ZETAF=ZETA 


HFL0914 




UEL=UL-G7XCU((ZETAF-ZETAL)/ZETAL 


HFL0915 




PPR=(ZETAF/ZETAR))^XNG 


HFL0916 




EVERR=PPR-1 . 


HFL0917 




SQRR=DSQRT(1 .+G6XEVERR) 


HFL0918 




USR=UR+CRGXEVERR/SQRR 


HFL0919 




DU=UEL-USR 


HFL0920 




IF (DABS(DU) .LE.UMIDA) GO TO 10 


HFL0921 




DUDZR = NGXCRGX(PPR/ZETAF)3((1 .+G9XEVERR)/SQRR3(X3 


HFL0922 




ZETA=ZETAF+DU/(DUDZR-DUDZL) 


HFL0923 




GO TO 11 


HFL0924 


10 


CONTINUE 


HFL0925 




USTAR=(UEL+USR)/2. 


HFL0926 




PSTAR=PPRXPR 


HFL0927 




CSTARL=CL + (UL‘*USTAR)/G7 


HFL0928 




RSTARL=GAMAXPSTAR/CSTARLXX2 


HFL0929 




GSTARL=CSTARLXRSTARL 


HFL0930 


C EQU. NO. 69.01 OF THE BOOK BY COURANT-FRI EDRICHS . 


HFL0931 




NWR = G11X(USTAR“UR)3(R0R 


HFL0932 




WR = WWR+DSQRT(GR5^X2+WWR3e)^2) 


HFL0933 




RSTARR=RORXWR/(WR-RORX(USTAR“UR) ) 


HFL0934 




GSTARR=DSQRT(GAMAXPSTARXRSTARR) 


HFL0935 




CSTARR=GSTARR/RSTARR 


HFL0936 
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ILE: HFL 



SCRIPT A1 



WRE=WR/ROR+UR 

UW(1)=UL-CL 

UW(2)=USTAR“CSTARL 

UW(3)=USTAR 

UW(A)=WRE 

UW(5)=WRE 

GO TO 5 

C THE CASE SS 

2 ITYPE=NCASE 
HELEML^.TRUE. 

HELEMR=,TRUE. 

21 N=N+1 

IF (N.GT.NMAX) GO TO 7003 

ZETAF=ZETA 

PF = ZETAF^)^NG 

PPL=PF/PL 

PPR=PF/PR 

EVERL=PPL-1 . 

EVERR=PPR-1 . 

SQRL = DSQRT(1 .+G65(EVERL) 

SQRR = DSQRT(1 .+G65(EVERR) 

USL = UL-CLG5^EVERL/SQRL 
USR=UR+CRG^EVERR/SQRR 
DU=USL-USR 

IF (DABS(DU) .LE.UMIDA) GO TO 20 

DUDZL = -NG5^CLG^(PPL/ZETAF)5((1 .+G95(EVERL )/SQRL^^3 

DUDZR= NG5(CRG5((PPR/ZETAF)5((1 .+G93^EVERR)/SQRR)^^3 

ZETA=ZETAF+DU/(DUDZR-DUDZL) 

GO TO 21 
20 CONTINUE 

USTAR=(USL+USR)/2. 

PSTAR = (PPU^PL+PPR5(PR)/2. 

WWR = G115((USTAR-UR)5(R0R 
WR=WWR+DSQRTCGR^5(2+WWR5()(2) 
WWL=-’Gllx(USTAR-UL)5(ROL 
NL = WWL + DSQRT(GL^5^2+WWL5(^2) 

RSTARL=ROL5(WL/(WL + ROL5((USTAR-UL) ) 

RSTARR = ROR^WR/(WR-ROR5((USTAR-UR) ) 

GSTARL = DSQRT(GAMA5(PSTAR^RSTARL) 

GSTARR = DSQRT(GAMA5(PSTAR5(RSTARR) 

CSTARL=GSTARL/RSTARL 

CSTARR=GSTARR/RSTARR 

WLE=-WL/ROL+UL 

WRE=WR/ROR+UR 

UW(1)=HLE 

UW(2)=WLE 

UW(3)=USTAR 

UWCA)=WRE 

UW(5)=WRE 

GO TO 5 

C THE CASE SE 

3 ITYPE=NCASE 
HELEML=.TRUE. 

HELEMR=. FALSE. 

31 N=N+1 

IF (N.GT.NMAX) GO TO 7003 
ZETAF=ZETA 

UER = UR+G75^CR)((ZETAF-ZETAR)/ZETAR 
PPL = (ZETAF/ZETAL)5(5(NG 
EVERL=PPL-1. 

SQRL = DSQRT(1 .+G6 5^EVERL) 

USL=UL-CLG^EVERL/SQRL 

DU=USL-UER 

IF (DABS(DU) .LE.UMIDA) GO TO 30 

DUDZL=-NG^CLG^(PPL/ZETAF)5((1 . +G95^EVERL )/SQRL5()^3 
ZETA=ZETAF+DU/(DUDZR-DUDZL) 

GO TO 31 
30 CONTINUE 

USTAR=(USL+UER)/2. 

PSTAR=PPL^PL 
CSTARR=CR-(UR‘USTAR)/G7 
RSTARR = GAMA^PSTAR/CSTARR^5(2 
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HFL 



SCRIPT A1 



GSTARR=CSTARRXRSTARR 
WWL=-G11^((USTAR-UL)XR0L 
WL=WNL + DSQRT.(GLX3C2+WWL^X2) 

WLE=-NL/R0L+UL 

RSTARL=ROLXWL/(WL+ROL^(USTAR-UL) ) 

GSTARL = DSQRT(GAMA)(PSTARXRSTARL) 

CSTARL=GSTARL/RSTARL 

UW(1)=WLE 

UW(2)=WLE 

UW(3)=USTAR 

UW(^)=USTAR+CSTARR 

UN(5)=UR+CR 

GO TO 5 

C THE CASE EE 
A ITYPE=NCASE 

HELEML=. FALSE. 

HELEMR=. FALSE. 

PSTAR = ZETA)(XNG 

USTAR = UL-G73(CU^(ZETA-ZETAL)/ZETAL 

CSTARL=CL+(UL-USTAR)/G7 

CSTARR=CR-(UR-USTAR)/G7 

RSTARL=GAMAXPSTAR/CSTARLXX2 

RSTARR = GAMA5(PSTAR/CSTARR)^^2 

GSTARL=RSTARLXCSTARL 

GSTARR=RSTARRXCSTARR 

UN(1)=UL-CL 

UW(2)=USTAR-CSTARL 

UW(3)=USTAR 

UW(4)=USTAR+CSTARR 

UW(5)=UR+CR 

N = 1 

GO TO 5 

5 CONTINUE 
DO 6 K=l,6 
NFLUX=K 

IF (UW(K) .GE.O. ) GO TO 61 

6 CONTINUE 
NFLUX=6 

61 CONTINUE 

MIN(I)=MIN(I) .OR.NCASE 
MIN(I)=MIN(I) .0R.SHIFT(N,3) 

MIN(I)=MIN(I) .0R.SHIFT(NFLUX,9) 

NC14(NCASE)=NC14(NCASE)+1 
CASEAV(NCASE)=CASEAV(NCASE)+N 
NF16(NFLUX)=NF16(NFLUX)+1 
IF(A(3) .NE.O. )GO TO 666 
IF(I.NE.2) GO TO 666 

PRINT 667, I.NFLUX, NCASE, PL , UL , ROL , PR, UR, ROR, USTAR, PSTAR, RSTARL , 

1 RSTARR, (KK,UW(KK),KK=1,6) 

667 FORMAT (/IX, ’I,NFLUX, NCASE= ' , 3I5/1X, ’ PL , UL , ROL , PR, UR, ROR= * , 6D12 . 4/ 

1 IX, ’USTAR,PSTAR,RSTARL,RSTARR=',4D13.4/ 

2 IX, *KK,UW(KK)=* ,6(I4,2X,D13.4)/) 

666 CONTINUE 

RETURN 

7001 CONTINUE 

PRINT 7101, PL,UL,PR,UR,ZETAL,ZETAR,SL,SR,NL,NR, I 
FORMATC//1X, 'FROM RIEMAN. AN IMPOSSIBLE CASE OF EXPANSION/SHOCK' 

1 //IX, *PL,UL,PR,UR=' ,4D25.14// 

2 IX, 'ZETAL,ZETAR,SL,SR=',4D25.14// 

3 IX, 'NL,NR,I=',3I10//) 

CALL S0F(7001) 

7002 CONTINUE 

PRINT 7102, ZETA,DUDZL,DUDZR,ZETAL,ZETAR,PL,UL,PR,UR,N,NCASE,I 
7102 F0RMATC//1X, 'FROM RIEMAN. NEGATIVE PRESSURE AT THE INTERSECTION' 

1 IX, 'OF L AND R EXPANSION BRANCHES'// 

2 IX, 'IT MEANS THAT A CAVITATION TENDS TO FORM. THIS', 

3 IX, 'POSSIBILITY IS EXCLUDED IN PRESENT VERSION'// 

4 IX, 'ZETA,DUDZL, DUDZR,ZETAL,ZETAR,PL,UL,PR,UR=' ,9D10.3// 

5 IX, 'N,NCASE,I=',3I10//) 

CALL SOFC7002) 

7003 CONTINUE 

PRINT 7103, I,N,NCASE,DU,UMIDA,EPS,PL,UL,PR,UR, 
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ILE: HFL 



SCRIPT A1 



1 ZETA,ZETAF,ZETAL,ZETAR, DUDZL,DUDZR HFL 1081 

7103 FORMATC//1X, ’FROM RIEMAN. NUMBER OF ITERATIONS EXCEEDED.’// HFL1082 

1 .IX, ’I,N,NCASE, DU,UMIDA, EPS=’ ,3I6,3D18 .6// HFL 108 3 

2 ■ IX, ’PL,UL,PR,UR,ZETA,ZETAF=’,6D18.10// HFL108A 

3 IX, ’ZETAL,ZETAR,DUDZL,DUDZR=’ , AD18 .10//) HFL 108 5 

CALL S0FC7003) HFL1086 

RETURN HFL1087 

ENH HFl 1 088 



SUBROUTINE MAGA ( L , I , MIN ) 

IMPLICIT REALX8(A-H,0“Z,$) 

DIMENSION MIN(L) 

COMMON /GAM/ GAMA , NG, MU2 , G1 , G2 , G3 , GA , G5 , G6 , G7 , G8 , G9 , G1 0 , G1 1 
,G12,G13,G1A,G15,G16,G17,G18,G19,G20 

REALX8 NG 
REALX8 MU2 

COMMON /STEPO/UL,PL,ROL,GL,UR,PR,ROR,GR,USTAR,PSTAR, 

1 RSTARL, RSTARR,GSTARL,GSTARR,WL,WR, 

2 CL , CR, CSTARL , CSTARR, UW( 6 ) 

3 ,ZL,ZR,ZSTARL,ZSTARR,NFLUX,HELEML,HELEMR 
LOGICAL HELEML,HELEMR 

COMMON /STEPl/DUIDT, DPIDT,DGIDTL,DGIDTR, DRIDTL , DRIDTR 

2 , ASTARL,ASTARR 

3 ,RAT,SH 

COMMON /GRADS/DUDXIL, DPDXIL,DGDXIL,DRDXIL, 

1 DUDXIR,DPDXIR,DGDXIR,DRDXIR 

COMMON /AB/AC50) 

REAL LU,LP,LRO 
DATA EPS/1. D-6/ 



1 



HFL1089 

HFL1090 

HFL1091 

HFL1092 

HFL1093 

HFL109A 

HFL1095 

HFL1096 

HFL1097 

HFL1098 

HFL1099 

HFLllOO 

HFLllOl 

HFL1102 

HFL1103 

HFLllOA 

HFL1105 

HFL1106 

HFL1107 

HFL1108 



C WE HERE SOLVE FOR THE TIME-DERI VATIVES ALONG THE CONTACT SURFACE, HFLlllO 
C NAMELY DUIDT,DPIDT. FROM THESE WE ALSO OBTAIN THE OTHER HFLllll 
C TIME-DERIVATIVES (SEE COMMON /STEPl/). HFL1112 
C WE COMPUTE THE COEFFICIENTS FOR TWO EQUATIONS FOR DUIDT,DPIDT. THESEHFL1113 
C ARE AAL3(DUIDT + BBLXDPIDT=CCL HFLlllA 
C AARXDUIDT+BBRXDPIDT=CCR HFL1115 



: if(sh.le.eps)rat=o. 

IF ( .NOT.HELEML) GO TO 12 

11 CONTINUE 
DP=PSTAR-PL 

IF (DABS(DP) .LT.EPSxpSTAR) GO TO 12 

DU=USTAR-UL 

Zl = l ./DP 

Z2=0 . 5/(PSTAR+MU2xPL) 

ZA=Z1-Z2 

ZB=Z1+MU2xZ2 

LU=-ROLXDUX(GAMA^PLXZB+0 .5)+WL 
LRO=0.5XDP/ROL 
LP=1 .+ZBXDP 
AA I = 1 +7A ¥nP 

3BL=-ZAxDU+WL/GSTARL 5(X2 
CCL = -LU^DUDXIL-LROB(DRDXIL-LPxDPDXIL 
CCL=CCL-WL^USTARXRAT/RSTARL 
1 +UL5(RATXDUx(GAMAxPLXZB+0 .5) 

GO TO 10 

12 CONTINUE 
A1=DUDXIL+DPDXIL/GL 
BETA=GSTARL/GL 

ASTARL=A1 + (G3/GL)X(CLxDGDXIL-GAB(DPDXIL)X(BETAXXG5-1 , 
AAL=1 . 

BBL=1 ./GSTARL 

CCL=-GSTARLXASTARL/DSQRT(BETA) 

GEOM = RAT5((( GAMA-1 . )xUL + 2.xCL)^ 

I (BETA^XG13-1 . )/(R0Lx(GAMA-3 . ) ) 

1 -A.^RAT^CLX(BETAXXG1A-1 . )/ ( ROLX( 3 . )^GAMA-5 . ) ) 
ASTARL=ASTARL-GEOM 
EVER1= GSTARU^GEOM/DSQRTCBETA) 

EVER2 = -RATxUSTAR3(CSTARL 
CCL=CCL+EVER1+EVER2 
GO TO 10 
10 CONTINUE 

IF ( .NOT.HELEMR) GO TO 22 
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E: HFL 



SCRIPT A1 



21 CONTINUE 
DP=PSTAR-PR 

IF (DABS(DP)-.LT.EPSXPSTAR) GO TO 22 

DU=U5TAR-UR 

Zl=l ./DP 

Z2 = 0.5/(PSTAR+MU2)CPR) 

ZA=Z1-Z2 

ZB=Z1+MU2^Z2 

LU = -ROR)^DUX(GAMA^PRXZB+0.5)-WR 
LRO=0 .5^DP/R0R 
LP = 1.+ZB5(DP 
AAR = 1 .+ZA5(DP 
BBR=-ZAXDU-NR/GSTARRXX2 
CCR = -LUB(DUDXIR-LR0XDRDXIR-LPXDPDXIR 
CCR = CCR+WR5^USTARXRAT/RSTARR 
1 +UR)^RAT^DUX(GAMAXPR5(ZB + 0 .5) 

GO TO 20 

22 CONTINUE 
A1=DUDXIR“DPDXIR/GR 
BETA=GSTARR/GR 

ASTARR = A1-(G3/GR)X(CRXDGDXIR-G^)^DPDXIR)B((BETA^XG5-1 . ) 
AAR=1 . 

BBR=-1 ./GSTARR 

CCR=G5TARR^ASTARR/DSQRT(BETA) 

GE0M = RATX(-(GAMA-1 . ))CUR+2.^CR)5((BETAXXG13-1 . ) 

1 /(RORXCGAMA-3. )) 

2 -4 . )^RATXCRX( BETAXXG14-1 . )/( RORX( 3 . )^GAMA-5 . ) ) 
ASTARR=ASTARR+GEOM 
EVER1=GSTARR5^GE0M/DSQRT(BETA) 

EVER2=RAT^USTARXCSTARR 
CCR=CCR+EVER1+EVER2 
GO TO 20 
20 CONTINUE 

DET=AAL^BBR-AARB(BBL 
DUIDT=(CCLXBBR“CCR^BBL)/DET 
DPIDT = -(CCLXAAR-CCR)^AAL)/DET 
DRIDTL = DPIDT/CSTARLB(X2 
DRIDTR=DPIDT/CSTARR)^)^2 

DGIDTL=0.5XGSTARU^(DPIDT/PSTAR+DRIDTL/RSTARL) 

DGIDTR = 0 .5XGSTARR)^(DPIDT/PSTAR+DRIDTR/RSTARR) 

RETURN 
END 
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SUBROUTINE FLUXEC L , I , MIN) 

IMPLICIT REAL^8(A-H,0-Z,$) 

DIMENSION MIN(L) 

COMMON /AB/AC50) 

EQUIVALENCE ( DT , A ( ^ ) ) , ( NCYC, A ( 12 ) ) 

COMMON /GAM/GAMA. NG.MU2,G1,G2,G3,G^,G5,G6,G7,G8,G9,G10,G11 
1 ,G12,G13,G14,G15,G16,G17,G18,G19,G20 

REALX8 NG 
REALX8 MU2 

COMMON /GRADS/ DUDXI L , DPDXI L , DGDXI L , DRDXIL , 

1 DUDXI R. DPDXI R,DGDXIR,DRDXIR 

COMMON /STEPO/UL,PL,ROL,GL,UR.PR,ROR,GR,USTAR,PSTAR, 

1 RSTARL,RSTARR,GSTARL, GSTARR, WL,WR, 

2 CL,CR,CSTARL,CSTARR,UW(6) 

3 , ZL , ZR, ZSTARL , ZSTARR, NFLUX , HEL EML , HEL EMR 
LOGICAL HELEML,HELEMR 

COMMON /STEP1/DUIDT,DPIDT,DGIDTL,DGIDTR,DRIDTL,DRIDTR 

2 ,ASTARL,ASTARR 

3 ,RAT,SH 



HFL1195 

HFL1196 
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COMMON/DETO/QDET,TC,RATE,PCJDET,RCJDET,UCJDET,DCJDET,PODET,ROODET HFL121^ 
COMMON /FI/FIH1,FIH2,FIH3,UXN,PXN,GXN,R0XN HFL 1215 

1 ,GIH HFL1216 

2 ,FIH^,ZMDOTL,ZMDOTR HFL1217 

COMMON /PRNT/ROX,PX,UX,ZX HFL1218 

ENTRY FLUXE1(L,I,MIN) HFL1220 

DT2=DT/2. HFL1222 

GO TO (1,2, 3, 4, 5, 6), NFLUX HFL1223 

1 CONTINUE HFL1224 
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ILE: HFL 



SCRIPT A1 



DUDXIX=DUDXIL 

DRDXIX=DRDXIL 

DPDXIX-=JDPDXIL 

DGDXIX = 0 .5^GL)^(DPDXIL/PL + DRDXIL/R0L) 

DUDTX=-DPDXIL 
DR0DTX=-R0Lxx2xDUDXIL 
DPDTX=-GLxx2^DUDXIL 
DRODTX=DRODTX-RATxROL^UL 
DPDTX=DR0DTXXCLX^2 
DGDTX=G6XGL^DPDTX/PL 
UX = UL 
PX = PL 
ROX=ROL 
ZX = ZL 
GX = GL 
GO TO 9 
6 CONTINUE 

DUDXIX=DUDXIR 

DPDXIX=DPDXIR 

DRDXIX=DRDXIR 

DGDXIX=0 .5xGRx(DPDXIR/PR+DRDXIR/R0R) 

DUDTX=-DPDXIR 

DPDTX=-GRxx2xDUDXIR 

DR0DTX=-R0Rxx2^DUDXIR 

drodtx=drodtx-ratxrorxur 

DPDTX = DR0DTX)^CRxx2 
DGDTX=G6^GR^DPDTX/PR 
UX = UR 
PX = PR 
ROX=ROR 
ZX = ZR 
GX = GR 
GO TO 9 
2 CONTINUE 

BETA0=(MU2^(UL/CL+G7))xx(l ./MU2) 

A1=DUDXIL+DPDXIL/GL 

A0=Al+(G3/GL)x(CL^DGDXIL-GAxDPDXIL)^(BETA0^^G5-l . ) 

EVERl=-( ( GAMA-1) XUL+2XCL)^(BETA0^XG1 3-1 . )/( GAMA-3 . ) 

EVER2=A .^CLX(BETA0XXG1A-1 . )/ ( 3 . xGAMA-5 . ) 

EVER=(EVER1+EVER2)xRAT/R0L 

A0=(A0+EVER) 

DUDAX=A0 

dpdax=glxbetaoxao 

C0=MU2x(UL+G7^CL) 

UX = CO 

ROX=GLXBETAO/CO 

ZX=ZSTARL 

PX=ROXXC0^^2/GAMA 

GX=ROXXCO 

DRODAX=ROL^BETAOxx(l ./GA) 

1 x( (DRDXIL/ROL-DPDXIL/(GAMA^PL))XDSQRT(BETAO) 

2 +AO/CO) 

DGDAX=BETAO^DSQRT(BETAO)x(DGDXIL-GAXDPDXIL/CL) 

1 +GA^ROL^AO^BETAOxx( 1 ./GA) 

DPDAX=DPDAX+RATxUX^COXDSQRT(BETAO) 

GA1=1 ./GA+0 .5 

DRODAX=(DRDXIL-DPDXIL/(CL^CL) ) xBETAOxxGAl+DPDAX/(COxCO) 
DGDAX = 0.5^GAMA^(PXXDRODAX+ROX)^DPDAX)/GX 
DUDBX=-CL^BETA0XX(-1 ./GA)/GA 
DPDBX=PL^BETA0^^MU2/G6 
DRODBX=ROLXBETAO^X(-MU2)/GA 
DGDBX=GL 
GO TO 9 
5 CONTINUE 

BETA0=(MU2X(-UR/CR+G7))^^(1 ./MU2) 

A1=DUDXIR-DPDXIR/GR 

A0=A1+(G3/GR)x(-CR^DGDXIR+GA^DPDXIR)x(BETA0^xG5-1 . ) 
EVER1=(-(GAMA-1 . ) XUR+2^CR) BETA0XXG13-1 . )/(GAMA-3. ) 

EVER2 = -A .xCRx(BETA0^xg1A-1 . )/ ( 3 . )^GAMA-5 . ) 

EVER=(EVER1+EVER2)^RAT/R0R 

A0=(A0+EVER) 

DUDAX=A0 
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HFL 



SCRIPT A1 



DPDAX=-GRXBETA0)^A0 
C0 = MU25((-UR+G75^CR) 

DR0DAX = -R0RXBETA05(X(1 ./G4) 

1 x(C-DRDXIR/ROR+DPDXIR/(GAMA)^PR))XDSQRT(BETAO) 

2 +AO/CO) 

UX=-CO 

ROX=GR5^BETAO/CO 
ZX=ZSTARR 
PX = R0X5^C0XX2/GAMA 
Gx= ROX^^C 0 

DGDAX = BETAO)(DSQRT(BETAO)x(-DGDXIR+GA5(DPDXIR/CR) 

1 +GAXROR)^AO)^BETAOX)((1 ./GA) 

DGDAX=-DGDAX 

DPDAX=DPDAX-RAT5(UXXC05^DSQRT(BETAO) 

GAl=l./GA+0.5 

DRODAX=(DRDXIR-DPDXIR/(CR5^CR)) XBETA0)^XG41 + DPDAX/ ( COXCO ) 

DGDAX = 0.5XGAMAX(PX5(DR0DAX+R0XXDPDAX)/GX 

DUDBX=CRXBETAOXX(-l ./GA)/G4 

DPDBX=PRXBETA0XXMU2/G6 

DRODBX = ROR5^BETA05^)^(-MU2)/GA 

DGDBX=GR 

GO TO 9 

3 CONTINUE 
UX=USTAR 
PX=PSTAR 
ROX=RSTARL 
ZX=ZSTARL 
GX=GSTARL 

DUDXIX=-DPIDT/GSTARLxx2 

DPDXIX=-DUIDT 

DUDXIX=DUDXIX-RATXUSTAR/RSTARL 
IF ( .NOT.HELEML) GO TO 32 

31 CONTINUE 

DRDXIX=(RSTARL/WL)XX2X(3.XDUIDT+DPIDTX(1 .+3.X(WL/GSTARL))^X2)/WL 

1 +DUDXIL5^WU^((GL/WL)5(X2+3. )+3.5^DPDXIL 

2 +DRDXILX(WL/R0L)XX2) 
EVER1=ULXRSTARLXX25^RAT5(( (GL/WL)5(X2+1 . )/(ROL5(WL) 

EVER2 = 2 . 5^RSTARLXUSTARXRAT/WL 
DRDXIX=DRDXIX+EVER1+EVER2 

GO TO 33 

32 CONTINUE 
BETA=GSTARL/GL 
SQB=DSQRT(BETA) 

DRODA = ROL^BETA3(X(l ./GA)x(SQBx(DRDXIL/ROL-DPDXIL/(GAMA5^PL)) 

1 +ASTARL/CSTARL) 

DRDXIX=DR0DA/SQB + DPIDT/(GSTARL^CSTARLX5^2) 

DPDA = GSTARLX(ASTARL + RATXUSTARXCSTARL/(GL5^ SQB)) 

G41=l./G4+0.5 

DR0DA=(DRDXIL-DPDXIL/(CLXCD) XBETAxxG41+DPDA/( CSTARLXX2) 
DRDXIX= DRODA/SQB + DPIDT/CGSTARLXCSTARL 5^5(2) 

33 CONTINUE 

DGDXIX = 0.55^GXX(DPDXIX/PX+DRDXIX/R0X) 

DUDTX=DUIDT 

DPDTX=DPIDT 

DGDTX=G6XGSTARL5^DPIDT/PSTAR 
DR0DTX=DPIDT/CSTARLXX2 
GO TO 9 

4 CONTINUE 
UX=USTAR 
PX=PSTAR 
ROX=RSTARR 
ZX=ZSTARR 
GX=GSTARR 

DUDXIX=-DPIDT/GSTARRXX2 
DUDXIX=DUDXIX-RATXUSTAR/RSTARR 
DPDXIX=-DUIDT 
IF ( .NOT.HELEMR) GO TO 42 
41 CONTINUE 

DRDXIX=(RSTARR/HR)X)^2)^(3.5(DUIDT-DPIDT5^(1 .+3 .X(WR/GSTARR)X3(2)/WR 

1 -DUDXIRXWRX((GR/WR)xx2+3. )+3 .xdPDXIR 

2 +DRDXIRX(WR/R0R)X5^2) 
EVER1=UR5(RSTARR5^x2xRATX((GR/WR)xx2+1 . )/(ROR5^WR) 
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ILE: 



HFL SCRIPT A1 



EVER2 = 2 .XRSTARR)(USTAR)^RAT/WR HFL 156 9 

DRDXIX=DRDXIX-EVER1-EVER2 HFL 1370 

GO TO HFL1571 

A2 CONTINUE HFL1372 

BETA=GSTARR/GR HFL 157 5 

SQB = D5QRT(BETA) HFL 157 A 

DR0DA = -R0R)^BETAX^^(1 ./GA ))^(SQB)^(-DRDXIR/ROR+DPDXIR/(GAMA^PR) ) HFL 157 5 

1 +ASTARR/CSTARR) HFL1576 

DRDXIX = DR0DA/SQB-DPIDT/(GSTARR)^C5TARR)^)^2) HFL 157 7 

DPDA = -GSTARR^(ASTARR+RAT^USTAR)^CSTARR/(GRX SQB)) HFL 157 3 

GA1=1 ./GA+0.5 HFL1579 

DRODA = (DRDXIR-DPDXIR/(CR)^CR)) ^BETA^^GAl + DPDA/ (CSTARR)()^2 ) HFL 158 0 

DRDXIX= DR0DA/SQB-DPIDT/(GSTARR^CSTARR)^)^2) HFL 1581 

A5 CONTINUE HFL1582 

DGDXIX = 0.5^GX)^(DPDXIX/PX+DRDXIX/R0X) HFL 1385 

DUDTX=DUIDT HFL158A 

DPDTX=DPIDT HFL1385 

DGDTX = G6)(G5TARR^DPIDT/PSTAR HFL 1386 

DR0DTX = DPIDT/CSTARR)^)^2 HFL 1587 

GO TO 9 HFL1588 

9 CONTINUE HFL1589 

RETURN HFL1590 

ENTRY FLUXE2(L,I,MIN) HFL1592 



FI1=R0X^UX 

FI2=R0XXUXXX2+PX 

FI2=FI2-PX 

FI5 = UX^(G12xPX+0 . 5XR0X)^UX)(^2+ZX^R0XXQDET) 

FIA=ZX^ROX)^UX 

ROUOO=ROX^UX 

GO TO(10,20,30,A0,50,60), NFLUX 
10 CONTINUE 
60 CONTINUE 

DFDXI1=DRDXIX^^UX+R0XXDUDXIX 

DFDXI2 = DRDXIX^UX)^^^2+2.^^R0X^UX)(DUDXIX+DPDXIX 

DFDXI2=DFDXI2-DPDXIX 

DFDXI3 = DUDXIX)^(G12)^PX+0.5^ROX)^UX)^X2) 

1 +UXX(G12^DPDXIX+0 . 5^DRDXIX)^UX)(X2+R0X^^UX)^DUDXIX) 

DFIDT1 = DR0DTX)^UX+R0X)CDUDTX 
DFIDT2 = DR0DTX)^UXX^2+2. ^ROXXUX^^DUDTX+DPDTX 
DFIDT2=DFIDT2-DPDTX 

DFIDT5 = DUDTX)^(G12)^PX+0 .5^R0X)CUX^^^2) 

1 +UXX(G12xDPDTX+0.5)^DR0DTX^UX^X2+R0X)^UX^DUDTX) 

FIDOT1 = -ROUOO^^DFDXI1 + DFIDT1 
FIDOT2=-ROUOO^DFDXI2+DFIDT2 
FIDOT5=-ROUOO^DFDXI3+DFIDT5 
UXDOT = -ROUOO)^DUDXIX+DUDTX 
PXDOT = -ROUOO)^DPDXIX+DPDTX 
GXDOT=-ROUOO^DGDXIX+DGDTX 
ROXDOT=~ROUOOXDRDXIX+DRODTX 
FIH1 = FI1 + DT2)^FID0T1 
FIH2=FI2+DT2XFID0T2 
GIH=PX+DT2^PXD0T 
FIH5 = FI3+DT2)^FID0T3 
FIHA=FIA 
UXN = UX+DT^^UXDOT 
PXN = PX+DT^^PXDOT 
GXN=GX+DTxGXDOT 
ROXN = ROX+DT^^ROXDOT 
GO TO 90 
20 CONTINUE 

EVO=GL)^DSQRT(BETAO) 

201 CONTINUE 

DFIDA1 = DR0DAX)^UX+R0X^DUDAX 

DFIDA2 = DR0DAX^UXXX2+2 . )(ROX^UXXDUDAX+DPDAX 

DFIDA2=DFIDA2-DPDAX 

DFIDA5 = DUDAX)^(G12)(PX+0 .5^R0X)^UX)()^2) 

1 +UX>^(G12XDPDAX+0 .5^DR0DAX)^UX^^2+R0X^UXxDUDAX) 

FIDOT1=-EVO)(DFIDA1 
FIDOT2 = -EVO)CDFIDA2 
FID0T3=-EV0XDFIDA5 
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SCRIPT A1 



E: HFL 



FIH1=FI1+DT2^FID0T1 

FIH2=FI2+DT25(FID0T2 

FIH5=FI5+DT2XFID0T3 

FIHA=FIA 

GA=DGDAX 

IF(NFLUX.EQ.5)GA=-GA 
DROUA=UX^DRODAX+ROXXDUDAX 
BETAPR = 0.5)^DSQRT(BETA0)^(GA-DROUA) 
FIH2=FIH2“DPDBXXBETAPR)^DT2 
UXDOT = -EV05^DUDAX+BETAPR^DUDBX 
PXDOT=-EVOXDPDAX+BETAPR^DPDBX 
GIH=PX+DT2XPXD0T 
GXDOT=-EVO^DGDAX+BETAPR^DGDBX 
ROXDOT=-EVO^DRODAX+BETAPR^DRODBX 
UXN=UX+DTXUXDOT 
PXN = PX+DT)^PXDOT 
GXN=GX+DT)^GXDOT 
ROXN=ROX+DTXROXDOT 
GO TO 90 
50 CONTINUE 

EVO=-GR^DSQRT(BETAO) 

GO TO 201 
30 CONTINUE 
40 CONTINUE 
GO TO 60 
90 CONTINUE 
RETURN 
END 
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